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Three  Dimensional  Transient  Analysis  of  Microstrip  Circuits  in  Multilayered  Anisotropic  Media 


Under  the  sponsorship  of  the  ONR  contract  N00014-90-J-1002  we  have  published 
12  refereed  journal  and  conference  papers. 

The  finite  difference-time  domain  (FD-TD)  technique  is  applied  to  the  solution 
of  Maxwell’s  equations.  A  computer  program,  which  can  be  used  to  simulate  and  study 
numerous  electromagnetic  phenomena,  is  developed  and  implemented  on  an  IBM  386  com¬ 
patible  personal  computer.  The  FD-TD  technique  is  a  useful  tool  for  students  in  electro¬ 
magnetics.  The  technique  is  flexible  and  can  be  applied  to  many  basic  EM  scattering  and 
radiation  problems.  Because  field  solutions  are  found  as  a  function  of  time,  visualization 
of  the  propagation  of  the  EM  fields  is  possible.  The  FD-TD  technique  is  implemented  for  a 
two-dimensional  rectangular  grid  in  conjunction  with  a  second-order  absorbing  boundary 
condition.  Both  E-  and  H-field  polarizations  are  analyzed.  Finite  objects  consisting  of 
dielectric,  magnetic  and  conducting  materials,  and  perfectly  conducting  infinite  ground 
planes  are  modeled.  Plane  wave  and  line  current  sources  are  implemented.  In  addition 
to  the  capability  of  animating  the  propagation  of  the  EM  fields,  radiation  and  scattering 
patterns  can  be  generated. 

A  methodology  developed  to  handle  dispersive  materials  in  the  time  domain  is 
extended  to  model  the  dispersive  characteristics  of  the  impedance  boundary  condition 
used  for  a  thin  layer  coating  over  perfect  conductors.  The  impedance  boundeiry  condition 
is  first  approximated  as  a  rational  function  of  frequency.  This  rational  function  is  then 
transformed  to  a  time  domain  equation,  resulting  in  a  partial  differential  equation  in 
space  and  time.  Discretization  of  the  time  domain  model  to  efficiently  handle  the  thin 
layer  coating  is  presented  in  the  context  of  the  finite-difference  time-domain  (FD-TD) 
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technique.  The  methodology  is  verified  by  solving  a  one-dimensional  problem  using  the 
FD-TD  technique  and  comparing  with  the  analytical  results. 

To  understand  the  physical  meaning  of  rational  reflection  coefficients  in  inverse¬ 
scattering  theory  for  optical  waveguide  design,  we  studied  the  relationship  between  the 
poles  of  the  transverse  reflection  coefficient  and  the  modes  in  inhomogeneous  dielectrics. 
By  using  a  stratified-medium  formulation  we  showed  that  these  poles  of  the  spectral  re¬ 
flection  coefficient  satisfy  the  same  equation  as  the  guidzmce  condition  in  inhomogeneous 
waveguides.  Therefore,  in  terms  of  wave  numbers,  the  poles  are  the  same  as  the  discrete 
modes  in  the  waveguide.  The  radiation  modes  have  continuous  real  values  of  transverse 
wave  numbers  £ind  are  represented  by  the  braiich  cut  on  the  complex  plane.  Based  on 
these  results,  applications  of  the  Gel’fand-Levitan-Marchenlro  theory  to  optical  waveguide 
synthesis  with  the  rational  function  representation  of  the  transverse  reflection  coefficient 
are  discussed. 

The  coupled-wave  theory  is  generalized  to  analyze  the  diffraction  of  waves  by  chiral 
gratings  for  arbitrary  angles  of  incidence  and  polarizations.  Numeric2il  residts  for  the 
Stokes  parameters  of  diffracted  Floquet  modes  versus  the  thickness  of  chiral  gratings  with 
various  chiralities  are  calculated.  Both  horizontal  and  vertical  incidences  ue  considered  for 
illustration.  The  diffracted  waves  from  chiral  gratings  are  in  general  elliptlcally  polarized; 
and  in  some  particular  instances,  it  is  possible  for  chiral  gratings  to  convert  a  linearly 
polarized  incident  field  into  two  nearly  circularly  polarized  Floquet  modes  propagating  in 
different  directions. 

A  general  spectral  domain  formulation  to  the  problem  of  radiation  of  arbitrary 
distribution  of  sources  embedded  in  a  horizontally  stratified  zirbitrary  magnetized  linear 
plasma  is  developed.  The  fields  are  obtained  in  terms  of  electric  and  magnetic  type  dyadic 
Green’s  functions.  The  formulation  is  considerably  simplified  by  using  the  kDB  system 
of  coordinates  in  conjunction  with  the  Fourier  transform.  The  distributional  sing\ilar 
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behavior  of  the  various  dyadic  Green’s  functions  in  the  source  region  is  investigated  and 
taken  into  account  by  extracting  the  delta  function  singularities.  Finally,  the  fields  in  any 
arbitrary  layer  are  obtained  in  terms  of  appropriately  defined  global  upward  zind  downward 
reflection  and  transmission  matrices. 

We  investigated  a  method  for  the  calculation  of  the  current  distribution,  resistance, 
and  inductemce  matrices  for  a  system  of  coupled  superconducting  transmission  lines  having 
finite  rectangular  cross  section.  These  calculations  allow  accurate  characterization  of  both 
high-Tc  and  low- Tc  superconducting  strip  transmission  lines.  For  a  single  stripline  geome¬ 
try  with  finite  ground  planes,  the  current  distribution,  resistance,  inductance,  and  kinetic 
inductance  are  calculated  as  a  function  of  the  penetration  depth  for  various  film  thickness. 
These  calculations  are  then  used  to  determine  the  penetration  depth  for  Nb,NbN,  and 
Y Ba2Cu2,0'i -X  superconducting  thin  films  from  the  measured  temperature  dependence 
of  the  resonant  frequency  of  a  stripline  resonator.  The  calculations  are  also  used  to  convert 
measured  temperature  dependence  of  the  quality  factor  to  the  intrinsic  surface  resistance 
as  a  function  of  temperature  for  a  Nb  stripline  resonator. 

The  electromagnetic  radiation  from  a  VLSI  chip  package  and  heatsink  structure 
is  analysed  by  means  of  the  finite-difference  time-domain  (FD-TD)  method.  The  FD- 
TD  algorithm  implemented  incorporates  a  multi-zone  gridding  scheme  to  accommodate 
fine  grid  cells  in  the  vicinity  of  the  heatsink  and  package  cavity  and  sparse  gridding  in 
the  remainder  of  the  computational  domain.  The  issues  pertaining  to  the  effects  of  the 
heatsink  in  influencing  the  overall  radiating  capacity  of  the  configuration  are  addressed. 

AnaJyses  are  facilitated  by  using  simplified  heatsink  models  and  by  using  dipole  elements  as  ' _ 

sources  of  electromagnetic  energy  to  model  the  VLSI  chip.  The  potential  for  enhancement 
of  spurious  emissions  by  the  heatsink  structure  is  examined.  For  heatsinks  of  typical 
dimensions,  resonance  is  possible  within  the  low  gigahertz  frequency  range.  _ 

Because  the  effects  of  diffraction  during  proximity-print  x-ray  lithography  are  of 
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critical  importance,  a  number  of  previous  researchers  have  attempted  to  calculate  the 
diffraction  patterns  and  minimum  achievable  feature  sizes  as  a  function  of  wavelength  emd 
gap.  Work  to  date  has  assumed  that  scalar  diffraction  theory  is  applicable-as  calculated,  for 
example,  by  the  Rayleigh-Sommerfeldformulation-and  that  Kirchhoff  boundary  conditions 
Ccin  be  applied.  Kirchhoff  boundary  conditions  assume  that  the  fields  (amplitude  and 
phase)  are  constant  in  the  open  regions  between  absorbers,  and  a  different  constant  in 
regions  just  under  the  absorbers  (i.e.,  that  there  are  no  fringing  fields).  An  x-ray  absorber 
is,  however,  best  described  as  a  lossy  dielectric  that  is  tens  or  hundreds  of  wavelengths 
tall,  and  hence  Kirchhoff  boundary  conditions  are  unsuitable.  We  have  instead  used  two 
numerical  techniques  to  calculate  accurate  diffracted  fields  from  gold  absorbers  for  two 
cases:  a  30  nm-wide  line  at  A  =  4.5  nm,  and  a  100  nm-wide  line  at  A  =  1.3  nm.  We 
show  that  the  use  of  Kirchhoff  boundary  conditions  introduces  unphysically  high  spatial 
frequencies  into  the  diffracted  fields.  The  suppression  of  these  frequencies-which  occurs 
naturally  without  the  need  to  introduce  an  extended  source  or  broad  spectrum-improves 
exposure  latitude  for  mask  features  near  0.1  fim.  and  below. 

In  order  to  understand  the  physical  meaning  of  rational  reflection  coefficients  in 
one-dimensional  inverse  scattering  theory  for  optical  waveguide  design,  we  have  studied 
the  relation  between  the  poles  of  the  transverse  reflection  coefficient  and  the  modes  in 
inhomogeneous  dielectrics.  By  using  a  stratified  medium  model  it  is  shown  that  these 
poles  of  the  reflection  coefficient  have  a  one-to-one  correspondence  to  the  discrete  modes, 
which  are  the  guided  and  leaky  modes.  The  radiation  modes  have  continuous  real  values  of 
transverse  wave  numbers  and  are  not  represented  by  the  poles  of  the  reflection  coefficient. 
Based  on  these  results,  applications  of  the  Gel’fand-Levitan-Marchenko  theory  to  optical 
waveguide  synthesis  with  the  rational  function  representation  of  the  transverse  reflection 
coefficient  are  investigated. 

In  compact  modules  of  high  performance  computers,  signal  transmission  lines  be- 
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tween  integrated  circuit  chips  are  embedded  in  multilayered  dielectric  medium.  These 
signal  lines  are  usually  placed  in  different  layers  and  run  perpendicular  to  each  other.  The 
interaction  between  the  orthogonal  crossing  lines  and  the  signzd  line  affects  its  propagation 
characteristics  and  may  cause  considerable  signal  distortion. 

The  interaction  of  a  pair  of  crossing  lines  in  isotropic  medium  has  been  studied  using 
a  time-dommn  approach,  where  coupling  is  described  qualitatively.  This  method  becomes 
computationedly  expensive  when  the  number  of  crossing  lines  increases.  With  many  identi¬ 
cal  crossing  strips  uniformly  distributed  above  the  signal  line,  the  transmission  properties 
zire  characterized  by  stopbands  due  to  the  periodicity  of  the  structure.  Periodic  struc¬ 
ture  have  been  investigated  using  frequency-dommn  methods.  Periodically  nonuniform 
microstrip  lines  in  an  enclosure  are  analyzed  on  the  basis  of  a  numerical  field  calculation. 
A  technique  based  on  the  network- analytical  formulism  of  electromagnetic  fields  has  been 
used  to  analyze  striplines  and  finlines  with  periodic  stubs.  The  propagation  chuacteristics 
of  waves  along  a  periodic  array  of  parallel  signal  lines  in  a  multilayered  isotropic  struc¬ 
ture  in  the  presence  of  a  periodically  perforated  ground  plane  and  that  in  a  mesh-plane 
environment  have  been  studied.  More  recently,  the  effect  of  the  geometrical  properties  on 
the  propagation  characteristics  of  strip  lines  with  periodic  crossing  strips  embedded  in  a 
shielded  one-layer  isotropic  medium  have  been  investigated.  Both  open  and  closed  mul¬ 
tilayered  uniaxially  anisotropic  structures  are  considered.  A  full-wave  analysis  is  used  to 
study  the  propagation  characteristics  of  a  microstrip  line  in  the  presence  of  crossing  strips. 
The  signal  line  and  the  crossing  strips  are  assumed  to  be  located  in  two  arbitrary  layers  of 
a  stratified  uniaxially  anisotropic  medium.  An  integral  equation  formulation  using  dyadic 
Green’s  functions  in  the  periodically  loaded  structure  is  derived.  Galerkin’s  method  is 
then  used  to  obtain  the  eigenvalue  equation  for  the  propagation  constant.  The  effects  of 
anisotropy  on  the  stopband  properties  are  investigated.  Numerical  results  for  open  and 
shielded  three-layer  uniaxially  anisotropic  media  are  presented. 
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For  microwave  integrated  circuit  applications,  the  characteristics  of  interconnects 
have  been  investigated  for  the  propagation  modes,  time  response,  crosstalk,  coupling, 
delay,  etc.  In  these  analyses,  it  is  assumed  that  quasi-TEM  modes  are  guided  along  the 
multiconductor  transmission  lines.  The  analysis  were  performed  for  arbitrary  number  of 
transmission  lines  where  the  load  and  the  source  conditions  were  presented  in  terms  of  the 
modal  reflection  and  transmission  coefflcient  matrices. 

To  perform  the  quasi-TEM  analysis,  the  capacitance  matrix  for  the  multiconductor 
transmission  line  has  to  be  obtained  first.  Both  the  spectral  and  the  spatial  domain 
methods  have  been  proposed  to  calculate  the  capaciteince  matrix.  In  the  spectral  domain 
methods,  two  side  walls  are  used  to  enclose  the  whole  transmission  line  structure,  and  the 
thickness  of  the  strip  lines  has  not  been  considered.  In  using  the  spatial  domain  method, 
the  structure  has  to  be  truncated  to  a  finite  extent  to  make  the  numerical  implementation 
feasible.  The  infinite  extent  of  the  structure  was  also  incorporated,  but  only  a  two-layer 
medium  was  considered. 

In  practical  microwave  integrated  circuits,  the  dielectric  loss  due  to  the  substrate 
and  the  conductor  loss  due  to  the  metallic  strips  are  also  studied  in  the  analysis  of  circuit 
performances. 

Based  on  the  scalar  Green’s  function,  a  set  of  coupled  integral  equations  is  obtained 
for  the  charge  distribution  on  the  strip  surfaces.  Pulse  basis  functions  emd  a  point-matching 
scheme  is  used  to  solve  numerically  the  set  of  integral  equations  for  the  charge  distribution, 
and  hence  the  capacitance  matrix.  The  duality  between  the  electrostatic  formulation  and 
the  magnetostatic  one  is  applied  to  calculate  the  inductance  matrix.  The  conductance 
matrix  is  obtained  by  using  the  duality  between  the  electrostatic  problem  and  the  current 
field  problem.  A  perturbation  method  is  used  to  calculate  the  resistance  matrix. 
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Finally,  a  transmission  line  analysis  is  derived  to  obtain  the  transfer  matrix  for  multi¬ 
conductor  uniform  lines,  which  significantly  reduces  the  effort  in  treating  the  load  and  the 
source  conditions.  Transient  responses  are  obtained  by  using  the  Fourier  transform.  The 
results  for  two  coupled  lines  are  obtained. 

With  the  ever  increasing  speed  and  density  of  modern  integrated  circmts,  the  need 
for  electromagnetic  wave  analysis  of  phenomena  such  as  the  propagation  of  transient  sig¬ 
nals,  especially  the  distortion  of  signal  pulses,  becomes  crucial.  One  of  the  most  important 
causes  of  pulse  distortion  is  the  frequency  dependence  of  conductor  loss,  which  is  caused 
by  the  “skin  effect”,  and  which  can  be  incorporated  into  the  circuit  models  for  trainsmis- 
sion  lines  as  frequency-dependent  resistance  and  inductance  per  unit  length.  Efficient  and 
accurate  algorithms  for  calculating  these  parameters  are  increasingly  important. 

We  have  developed  a  hybrid  cross-section  finite  element/coupled  integral  equation 
method.  The  technique  is  a  combination  of  a  cross-section  finite  element  method,  which  is 
best  for  high  frequencies.  An  interpolation  between  the  results  of  these  two  methods  gives 
very  good  results  over  the  entire  frequency  range,  even  when  few  basis  functions  are  used. 

In  the  cross-section  method,  we  divide  each  conductor  into  triangular  patches  and 
choose  one  of  the  patches  from  the  return  conductor  to  be  our  reference.  We  then  calculate 
the  resistance  and  inductance  matrices  for  the  patches.  Using  two  conditions  on  the  system, 
that  the  total  current  in  each  wire  is  the  sum  of  the  currents  in  the  patches,  and  that  the 
voltage  on  each  patch  in  a  wire  must  be  the  same  (no  transverse  currents),  we  can  reduce 
the  matrices  for  the  patches  to  the  matrices  for  the  wires.  In  the  Weeks  method,  the 
patches  are  rectangles,  and  the  quadruple  integral  is  done  quite  easily  in  closed  form. 
However,  it  is  also  possible  to  evaluate  the  quadruple  integral  in  closed  form  for  triangular 
patches,  although  the  mathematics  leading  to  this  result  is  quite  involved,  and  the  final 
form  of  the  answer  is  complicated.  We  therefore  use  triangular  patches  as  the  most  flexible 
means  of  modelling  conductors  with  arbitrary  cross-sections;  polygons  are  covered  exactly. 
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and  we  are  able  to  model  quite  closely  other  shapes,  such  as  circles. 

As  frequency  increases,  the  need  to  keep  the  uniform  current  approximation  valid  in 
the  patches  requires  either  the  addition  of  many  more  patches  as  the  skin  depth  decreases, 
or  a  redistribution  of  the  existing  patches  to  the  surface,  where  the  current  is.  However, 
changing  the  distribution  of  patches  makes  it  necessary  to  recalculate  the  resistance  zmd 
inductance  matrices  of  the  patches,  thus  increasing  the  computation  time.  Since  we  use  a 
surface  integral  equation  method  for  high  frequencies,  we  do  not  change  the  distribution 
of  the  triajigular  patches  for  the  cross-section  method  as  we  increase  the  frequency. 

For  high  frequencies,  we  use  a  coupled  surface  integral  equation  technique.  Under 
the  quasi-TEM  assumption,  the  frequency-dependent  resistance  and  inductance  result  from 
the  power  dissipation  and  magnetic  stored  energy,  which  can  be  calculated  by  solving  a 
magnetoquasistatic  problem,  with  the  vector  potential  satisfying  Laplace’s  equation  in  the 
region  outside  all  the  conductors.  The  resistance  and  inductance  are  usually  given  by 
integreds  of  these  field  quantities  over  the  cross-sections  of  the  wires,  but  by  using  some 
vector  identities  it  is  possible  to  convert  these  expressions  to  integrals  only  over  the  surfaces 
of  the  wires.  These  expressions  conteiin  only  the  current  at  the  surface  of  each  conductor, 
the  derivative  of  that  current  normal  to  the  surface,  and  constants  of  the  vector  potentijd. 
A  coupled  integral  equation  is  then  derived  to  relate  these  quantities  through  Laplace’s 
equation  and  its  Green’s  function  outside  the  conductors  and  the  diffusion  equation  and  its 
Green’s  function  inside  the  conductors.  The  method  of  moments  with  pulse  basis  functions 
is  used  to  solve  the  integral  equations.  This  method  differs  from  previous  work  in  that  the 
calculation  of  resistance  and  inductance  is  based  on  power  dissipation  and  stored  magnetic 
energy,  rather  than  on  impedance  ratios.  It  will  therefore  be  more  easily  extended  to 
structures  where  non-TEM  propagation  can  occur. 

For  the  intermediate  frequency  range,  where  the  conductors  are  on  the  order  of  the 
skin  depth,  were  found  it  very  efficient  to  interpolate  between  the  results  of  the  cross- 
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section  and  surface  methods.  The  interpolation  function  was  based  on  the  average  size  of 
the  conductors,  measured  in  skin  depths,  and  was  of  the  form  1/(1  +  0.16a^/^^),  where  it 
a  is  the  average  cross-section  of  the  conductors,  and  S  is  the  skin  depth. 
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ABSTRACT 

The  finite  difference-time  domain  (FD-TD)  technique  is  applied  to  the  solution  of  Maxwell's 
equations.  A  computer  program,  which  can  be  used  to  simulate  and  study  numerous 
electromagnetic  phenomena,  is  developed  and  implemented  on  an  IBM  386  compatible 
personal  computer.  The  FD-TD  technique  is  a  useful  tool  for  students  in  electromagnetics. 
The  technique  is  flexible  and  can  be  applied  to  many  basic  EM  scattering  and  radiation 
problems.  Because  field  solutions  are  found  as  a  function  of  time,  visualization  of  the 
propagation  of  the  EM  fields  is  possible.  The  FD-TD  technique  is  implemented  for  a  two- 
dimensional  rectangular  grid  in  conjunction  with  a  second-order  absorbing  boundary  con¬ 
dition.  Both  E-  and  H-field  polarizations  are  analyzed.  Finite  objects  consisting  of  dielectric, 
magnetic  and  conducting  materials,  and  perfectly  conducting  infinite  ground  planes  are 
modeled.  Plane  wave  and  line  current  sources  are  implemented.  In  addition  to  the  ca¬ 
pability  of  animating  the  propagation  of  the  EM  fields,  radiation  and  scattering  patterns 
can  be  generated.  ®  1992  john  Wiley  &  Sons.  Inc. 


1.  INTRODUCTION 

Tools  which  students  can  use  as  aids  for  under¬ 
standing  and  visualizing  electromagnetics  are  scarce. 
In  an  attempt  to  develop  such  a  tool,  the  finite  dif¬ 
ference-time  domain  (FD-TD)  technique  is  used  to 
solve  Maxwell’s  equations  [1-8].  This  method  is 
especially  suited  for  visualizing  electromagnetic 
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phenomena,  since  electric  and  magnetic  fields  are 
calculated  everywhere  within  a  computational  do¬ 
main  as  a  function  of  time.  Thus,  it  is  possible  to 
observe  how  electromagnetic  fields  propagate 
through  space  with  time.  By  facilitating  the  visual¬ 
ization  of  various  electromagnetic  phenomena,  the 
FD-TD  code  can  be  very  useful  in  developing  in¬ 
tuition  regarding  electromagnetics. 

The  first  step  in  applying  the  FD-TD  technique 
involves  approximating  Maxwell’s  equations  in  dif¬ 
ferential  form  by  center  differences  in  space  and 
time.  The  locations  at  which  electric  and  magnetic 
fields  are  calculated  are  positioned  on  some  sort  of 
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grid.  For  two-dimensional  problems,  two  typical 
grids  are  a  rectangular  [1]  and  a  triangular  mesh 
[4],  A  significant  advantage  of  using  rectangular 
grids  over  triangular  grids  is  greater  simplicity. 
However,  since  the  scattering  object  is  discretized 
on  a  rectangular  grid,  curves  and  slanted  lines  are 
approximated  by  staircases. 

To  begin  the  FD-TD  simulation,  all  the  fields 
within  the  computational  domain  are  initially  set 
to  zero.  At  each  time  step,  the  electric  fields  are  cal¬ 
culated  in  terms  of  the  electric  and  magnetic  fields 
of  the  previous  time  step  using  the  difference  equa¬ 
tions  obtained  earlier  from  Maxwell’s  equations. 
Next,  the  magnetic  fields  are  calculated  in  a  similar 
manner.  Boundary  conditions  are  enforced  at  the 
outer  boundary  of  the  computational  domain  and 
at  all  dielectric  and  conducting  interfaces.  At  the 
outer  boundary,  a  second-order  absorbing  boundary 
condition  is  utilized  in  order  to  simulate  unbounded 
space  beyond  the  computational  domain  [9].  The 
tangential  electric  fields  are  set  to  zero  at  the  con¬ 
ducting  surfaces,  and  the  tangential  electric  and 
magnetic  fields  are  kept  continuous  at  the  dielectric/ 
magnetic  boundaries.  The  implementation  of  the 
excitation  source  typically  requires  that  certain 
electric  and  magnetic  fields  be  updated  at  each  time 
step.  These  steps  are  essentially  repeated  until  steady 
state  is  reached  for  a  sinusoidal  excitation,  or  until 
all  the  transient  scattered  fields  have  propagated  out 
of  the  computational  domain  for  a  Gaussian  pulse 
excitation. 

In  order  to  keep  the  program  simple,  only  two- 
dimensional  problems  are  examined  using  a  uni¬ 
form  rectangular  grid.  Despite  the  fact  that  only  two- 
dimensional  geometries  are  considered,  valuable 
insight  can  be  gained  through  the  observation  of 
FD-TD  simulations  of  electromagnetic  phenomena. 
Also,  in  order  to  allow  wide  distribution  of  this  pro¬ 
gram,  the  program  is  written  for  IBM  386  compat¬ 
ible  personal  computers.  Both  electric  and  magnetic 
field  polarizations  are  implemented  and  examined. 
Scatterers  can  have  arbitrary  geometries  of  finite  size, 
which  consist  of  perfect  conductors,  dielectrics  with 
finite  conductivity,  and  magnetic  materials  with  fi¬ 
nite  magnetic  conductivity.  In  addition,  perfectly 
conducting  infinite  ground  planes  are  considered. 
The  excitation  sources  include  plane  waves,  and 
single  or  multiple  line  current  sources.  These  exci¬ 
tation  sources  can  be  either  sinusoidal  or  Gaussian 
pulses  in  time.  Since  field  values  are  calculated  ev¬ 
erywhere  within  the  computational  domain,  they 
can  be  stored  and  displayed  using  color  plots.  By 
displaying  the  fields  for  sequential  time  steps  in  rapid 
succession,  the  propagation  of  the  fields  can  be  an¬ 


imated.  By  storing  tangential  electric  and  magnetic  f 

field  values  over  a  closed  surface  enclosing  the  scat- 
terer.  it  is  possible  to  calculate  radiation  patterns 
using  Huygens’  principle. 


2.1.  Maxwell's  Equations  in  Rectangular 
Coordinates 

In  implementing  the  finite  difference-time  domain 
technique.  Maxwell’s  equations  must  be  discretized 
in  space  and  time.  In  this  case.  Yee’s  lattice,  which 
is  a  rectangular  grid,  is  utilized  because  of  its  sim¬ 
plicity  [  1  ] .  Maxwell’s  equations  in  vector  differential 
form  in  an  isotropic  and  homogeneous  dielectric 
and  magnetic  material  are: 

dB 

V  X  £  =  -  —  -  a,„H  ( I ) 

at 

-  dD 

^^n  =  —  +  aE  (2) 

^■D  =  p  (3) 

r-fi  =  o  (4) 


B  =  pfl  =  ( 5 ) 

D  =  lE  =  (rfoE  (6) 

where  pr  and  are  relative  permeability  and  per¬ 
mittivity,  and  <7  and  <7,„  are  electric  and  magnetic 
conductivities,  respectively.  Maxwell’s  divergence 
equations.(3)and(4), are  satisfied  in  the  finite  dif¬ 
ference  scheme  by  applying  appropriate  initial  and 
boundary  conditions.  The  actual  difference  equa¬ 
tions  used  in  the  FD-TD  technique  are  based  upon 
Maxwell’s  curl  equations,  ( 1 )  and  (2).  and  the  con¬ 
stitutive  relations,  (5)  and  (6). 

In  a  rectangular  coordinate  system.  Equations 
( 1 ),  (2),  (5),  and  (6)  can  be  written  as  the  following 
set  of  scalar  equations. 


dlE 

dEy 

BE. 

- 

(7) 

dz 

By 

dHy 

BE. 

BEy 

(8) 

Bx  ' 

'  Bz  ' 

-  0,„Hy 

dH. 

BE, 

BE, 

-  a,„H: 

(9) 

By  ■ 

Bx 

and  the  constitutive  relations  arc: 
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for  the  £-fieId  polarization. 

The  following  notation  for  any  function  of  time 
and  space  will  be  used  in  the  finite  difference  equa¬ 
tions. 


For  two-dimensional  problems,  which  are  as¬ 
sumed  to  be  uniform  and  infinite  in  the  y  direction, 
all  the  partial  derivatives  with  respect  to  y  are  equal 
to  zero.  Maxwell’s  equations  in  arbitrary  homoge¬ 
neous,  isotropic  media,  for  two-dimensional  prob¬ 
lems  in  a  rectangular  coordinate  system,  decouple 
into  the  //-field  polarization. 


dHy  dE,  dEx 
dt  dx  dz 
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dt  dx 
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and  the  £-field  polarization. 
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In  free  space,  where  (r  =  Hr=  1  and  a  = 
Equations  ( 13)-(  18)  simplify  to 

“  0, 
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for  the  //-field  polarization  and 

dEy  _  dHx  dH, 
dt  dz  dx 

(22) 
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The  partial  derivatives  in  space  and  time,  within 
Maxwell’s  equations,  are  approximated  using  center 
differences. 


dm 


m  +  AU2)-m-  AU2) 


(26) 


The  electric  and  magnetic  field  components  are  in¬ 
terlaced  in  time,  and  are  calculated  in  a  leap-frog 
manner  (i.e.,  first  the  electric  fields  are  calculated, 
then  the  magnetic  fields  are  calculated,  and  the  se¬ 
quence  is  repeated).  The  electric  and  magnetic  field 
components  are  interlaced  spatially  a  half-grid  cell 
apart. 

Only  the  //-field  polarization  will  be  treated  ex¬ 
plicitly  because  the  £-field  polarization  can  be 
treated  easily  through  the  use  of  duality. 


2.2.  Finite  Difference  Equations  for  the 
H-Field  Polarization 


2.2.1.  Treatment  of  Free  Space.  For  the //-field  po¬ 
larization,  the  finite  difference  equations  in  free 
space  are  obtained  by  applying  center  differences  to 
Equations  ( 19)-(2I ). 


X  £;^/  +  ^  2  j  ~  E:^i,  k  + 


(27) 
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ET'  ^  ,  /cj 

£r'(/,/c  +  ^)  =  £?(/,/c  +  l) 

5. ‘■  +  ^) 


where 


At  =  cA/  =  7==  (’* ) 

Vno<o 

The  locations  of  the  field  components  on  a  unit 
cell  of  the  rectangular  grid  for  the  //-field  polariza¬ 
tion  are  shown  in  Figure  1. 

2.2.2.  Treatment  of  Dielectric  and  Magnetic  Ma¬ 
terial.  The  finite  difference  equations  in  an  isotro¬ 
pic,  homogeneous,  dielectric,  and  magnetic  material 


0.  *)  Ex 

Figure  1  Field  Locations  on  Rectangular  Grid  for  H- 
Field  Polarization. 


are  obtained  by  discretizing  Equations  ( 13)-(  15) 
for  the  //-field  polarization.  The  derivatives  are  re¬ 
placed  by  center  differences,  and  the  electric-  and 
magnetic-field  terms  involving  electric  and  magnetic 
conductivities  are  approximated  by  using  an  average 
of  the  field  values  at  a  half  time  step  before  and  after 
the  desired  time. 
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2.2.3.  Treatment  of  Interfaces  Between  Two  Me¬ 
dia.  The  interface  between  any  two  media  always 
occurs  at  integer  nodes  (/,  k)  (Figure  2).  At  the 
interface  between  two  dielectric  and  magnetic  ma¬ 
terials,  the  finite  difference  equation  for  the  tangen¬ 
tial  electric  field  (i.e.,  £»)  must  be  treated  carefully 
(Figure  3).  Using  the  integral  form  of  Equation  (2) 
to  calculate  the  electric  field  tangential  to  the  inter¬ 
face,  ,  yields 


Figure  3  Interface  Between  Two  Dielectric  and  Magnetic 
Media  for  //-Field  Polarization. 


Changing  Equation  ( 35 )  to  a  finite  difference  equa¬ 
tion  yields 

ET'  ^  +  5  .  A:) 

^  '2(e,i  +  <,;)  -  (g|  -h  g2)>;oAT' 

[2(«,i  +  «,2)  +  (<r,  +  a2)vo^T, 


_ 4Ar/A: _ 

2(«,i  +  «,:)  -  (ffi  +  <r2)77oAr 


(36) 


2.2.4.  Treatment  of  Perfect  Electric  Conduc¬ 
tor.  The  boundary  condition  at  perfect  electric 
conductors  states  that  the  tangential  electric  field 
must  be  zero. 

nX£  =  0  (37) 

The  field  locations  for  the  //-field  polarization  with 
respect  to  the  integer  nodes  (/,  k)  are  shown  in  Fig¬ 
ure  1.  From  Figures  1  and  2,  it  is  clear  that  the 
electric  fields,  which  are  calculated  at  points  on  the 
surface  of  the  perfect  electric  conductor,  are  always 
tangential  to  the  surface.  Thus,  in  the  finite  differ¬ 
ence-time  domain  scheme,  the  boundary  condition 
at  the  perfect  electric  conductor  can  be  satisfied  by 
simply  setting  these  electric  fields  equal  to  zero  at 
each  time  step. 


•  (i-l.i+l)  •(i,it+l)  •  O+l.t+I) 


Figure  2  Positioning  of  a  Media  on  the  Rectangular 
Grid. 


2.2.5.  Treatment  of  Perfect  Magnetic  Conduc¬ 
tor.  At  a  perfect  magnetic  conductor,  the  tangential 
magnetic  field  must  be  zero,  i.e., 

nXH  =  0  (38) 

This  boundary  condition  is  satisfied  when  the  tan¬ 
gential  electric  field  at  the  surface  of  the  perfect 
magnetic  conductor  is  calculated.  One  way  of  ob¬ 
taining  a  difference  equation  for  the  tangential  elec¬ 
tric  field  ( i.e.,  E^)  is  to  simplify  the  original  problem 
(Figure  4)  by  using  image  theory  (Figure  5).  The 
difference  equation  for  £j  at  the  interface  between 
a  perfect  magnetic  conductor  and  a  dielectric  and 
magnetic  material  is  given  by 
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Figure  4  Interface  Between  a  Media  and  a  Perfect  Mag¬ 
netic  Conductor. 


if 


tr  -*■  fir 

a  -*•  a,„ 

<r,„  -*■  <7 

All  Other  variables  remain  the  same. 


(45) 

(46) 

(47) 

(48) 


2.4.  Stability  and  Accuracy 

The  choices  of  Ajc,  Az,  and  At,  are  motivated  by 
reasons  of  accuracy  and  stability  [3,  6].  In  general, 
to  obtain  accurate  results.  Ax  and  Az  must  be  a 
small  fraction  (~io)  of  the  smallest  wavelength  in 
any  media  expected  in  the  model  or  of  the  smallest 
dimension  of  the  scatterer.  For  this  problem,  the 
spacing  in  the  x  and  z  directions  will  be  equal  in 
size. 
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2Ar 
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(39) 


2f/r''-  X 


A.v  =  Az  =  A  (49) 


To  ensure  the  stability  of  this  time-stepping  algo¬ 
rithm.  At  must  satisfy  the  following  condition  [3]. 


At  £ 


(50) 


2.5.  Implementation  of  Excitation 
Source 


The  implementation  of  perfect  magnetic  con¬ 
ductors  is  of  particular  importance.  By  duality,  the 
perfect  magnetic  conductor  for  the  //-field  polar¬ 
ization  becomes  the  perfect  electric  conductor  for 
the  £-field  polarization. 

2.3.  Finite  Difference  Equations  for  the 
E-Field  Polarization 

All  the  difference  equations  needed  for  the  £-fieid 
polarization  can  be  easily  obtained  through  the  use 
of  duality.  The  difference  equations  used  for  the  //- 
field  polarization  can  be  directly  applied  for  the  £- 
field  polarization  with  the  following  substitutions. 


£—  /? 

(40) 

//—  -£ 

(41) 

fio  -*•  fo 

(42) 

<0  “*■  fiO 

(43) 

Vo  -*•  Uvo 

(44) 

The  two  basic  types  of  sources,  which  are  plane 
waves  and  line  current  sources,  are  treated  differ- 


Figure  5  Image  Problem  of  Interface  Between  a  Media 
and  a  Perfect  Magnetic  Conductor. 
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Figure  6  Absorbing  Boundary,  Computational  Domain 
and  Closed  Surface  on  Which  Fields  are  Sampled  With  a 
Finite  Scatterer. 


ently.  When  the  excitation  source  is  assumed  to  be 
either  a  sinusoidal  or  Gaussian  pulse  plane  wave, 
the  computational  domain  is  separated  into  two  re¬ 
gions  in  order  to  facilitate  the  treatment  of  the  ex¬ 
citation  source  (Figures  6  and  7). 

For  a  finite  scatterer,  within  the  inner  region,  total 
fields  are  calculated,  while  in  the  outer  region,  only 
scattered  fields  are  calculated.  The  scattered  fields 
are  defined  to  be  the  difference  between  the  total 
fields  and  the  incident  plane  wave, 

HAx,  z,  /)  =  z,l)-  Hi(x,  z,l)  (51) 

E,(x,  z,  t)  =  E,{x,  z,  t)  -  £,(x,  2,  t)  (52) 

The  incident  fields  are  subtracted  from  field 
quantities  just  within  the  inner  regions,  when  they 
are  used  to  calculate  field  quantities  just  beyond  the 


inner  region.  The  incident  fields  arc  added  to  field 
quantities  just  outside  the  inner  regions,  when  they 
are  used  to  calculate  field  quantities  just  within  the 
inner  region  [10], 

For  the  //-field  polarization,  some  sample  finite 
difference  equations  in  the  vicinity  of  the  interface 
between  the  two  regions  are  (Figure  8): 

(.4,  A:  4) 

+  (53) 

>?oA-VL  \  2) 

-  k  +  ij 

-[£?,(/ +  1.  a) a)]) 

£"■(/ 4.  *).£,-(/ A-) 

(54) 
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Figure  7  Absorbing  Boundary,  Computational  Domain 
and  Surface  on  Which  Fields  are  Sampled  With  an  Infinite 
Ground  Plane  Geometry. 


Outer  Region: 
Scattered  Fields 


Inner  Region: 
Total  Fields 


Figure  8  Interface  Between  Inner  and  Outer  Regions. 
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(55) 


For  an  infinite  ground  plane  configuration,  the 
computational  domain  is  basically  divided  into  the 
region  above  the  ground  plane  and  the  region  below 
the  ground  plane.  For  plane-wave  illumination  from 
the  upper  half  space,  total  fields  are  calculated  in 
the  lower  region,  while  in  the  upper  region,  only 
scattered  fields  are  calculated.  In  this  case,  the  scat¬ 
tered  fields  are  defined  to  be  the  total  fields  minus 
the  incident  plane  wave  and  the  plane  wave  that 
would  have  been  reflected  from  a  uniform  infinite 
ground  plane. 


//d.v,  r,  /)  =  H,(x,  t) 

-  /)  -  HAx, 


(56) 

/) 


£,(.v,  /)  =  £,(.v,  t) 

-  t)  -  £,(.v,  r,  I) 


(57) 


The  removal  of  the  reflected  plane  wave  from  the 
upper  half  space  is  essential,  because  the  reflected 
plane  wave  is  infinite  in  size  and  could  not  otherwise 
be  adequately  modeled  within  the  finite  computa¬ 
tional  domain.  For  plane-wave  illumination  from 
the  lower  half  space  the  total  and  scattered  fields  are 
simply  calculated  in  the  opposite  regions.  The  im¬ 
plementation  of  a  plane  wave  is  basically  identical 
to  the  implementation  for  a  finite  scattered  except 
with  the  slightly  different  definition  of  the  scattered 
field. 

The  implementation  of  line  current  sources  is 
straightfonvard.  In  this  situation,  total  fields  are  cal¬ 
culated  within  the  entire  computational  domain. 
Line  current  sources  are  essentially  point  sources  in 
this  two-dimensional  simulation.  For  the  £-field 
polarization,  an  electric  line  current  in  the  r  direc¬ 
tion  is  the  source,  while  for  the  //-field  polarization 
the  source  is  a  magnetic  line  current  in  the  t' direc¬ 
tion.  For  a  sinusoidal  excitation,  the  value  of  each 
line  current  source  is  added  to  the  value  of  either 
£,  or  Hy.  depending  on  the  polarization,  at  a  single 
point  at  every  time  step.  Alternatively,  for  a  Gauss¬ 
ian  pulse  excitation,  the  value  of  the  appropriate 
field  component  is  set  to  be  equal  to  the  line  current 
source.  The  advantage  of  adding  field  values  instead 
of  imposing  field  values  is  that  fields  can  propagate 


through  the  current  source.  The  advantage  of  im¬ 
posing  field  values  is  that  the  value  of  the  field  at 
that  point  is  known  for  all  time. 

The  actual  amplitude  of  the  current  is  not  strictly 
determined  by  the  amplitude  of  the  adjustment  to 
a  field  quantity  at  a  particular  node.  The  fields  near 
a  point  source  are  singular.  Thus,  in  this  FD-TD 
scheme,  the  adjustment  of  a  field  value  at  a  single 
point  does  not  correspond  to  a  two-dimensional 
point  source  (i.e.,  a  line  current  source  with  zero 
cross  section ),  rather  this  implementation  more  ac¬ 
curately  models  a  line  current  source  with  a  finite 
nonzero  cross  section.  One  way  of  obtaining  the 
amplitude  of  the  line  current  source  is  by  performing 
some  sort  of  calibration  [8].  In  particular,  the  pre¬ 
diction  obtained  using  the  FD-TD  scheme  with  a 
particular  grid  spacing,  for  the  power  radiated  by  a 
line  current  source  in  free  space,  can  be  compared 
to  the  closed-form  solution.  The  power  radiated  in 
the  FD-TD  scheme  can  be  obtained  by  integrating 
the  time-average  Poynting  power  density  over  some 
closed  surface.  The  time-average  Poynting  power 
density  can  be  found  from  the  electric  and  magnetic 
fields  produced  by  the  FD-TD  simulation.  The  am¬ 
plitude  of  the  field  quantity  added  at  a  single  point 
for  a  particular  grid  spacing  can  be  related  to  the 
amplitude  of  the  current  in  the  line  source.  This 
calibration  factor  will,  in  general,  be  a  function  of 
the  grid  spacing  [8],  This  calibration  factor  is  also 
applicable  to  line  current  sources  which  interact  with 
each  other  and/or  arbitrary  media,  so  long  as  the 
grid  spacing,  which  is  used,  is  the  same  as  the  spacing 
used  in  performing  the  calibration. 

The  presence  of  reactive  fields  generated  by  a  line 
current  source,  which  decay  rapidly  further  away 
from  the  source,  can  cause  inaccuracies  in  calculat¬ 
ing  field  values  in  the  vicinity  of  the  line  source  [7]. 
Hence,  interaction  between  a  line  current  source 
and  some  arbitrary  media  may  not  be  adequately 
modeled  when  the  media  is  close  enough  to  be  af¬ 
fected  by  the  reactive  fields  generated  by  the  line 
current  source.  One  simple  remedy  to  this  problem 
involves  using  a  finer  mesh  in  the  vicinity  of  the 
line  current  source  [7]. 

2.6.  Absorbing  Boundary  Condition 

An  absorbing  boundary  condition  (ABC)  at  the 
outer  boundary  is  needed  to  make  the  computa¬ 
tional  domain  finite,  and  to  simulate  unbounded 
space  beyond  the  computational  domain  ( Figures 
6  and  7)  [9-12].  The  absorbing  boundary  condition 
used  for  this  problem  is  the  second-order  approxi¬ 
mation  derived  by  Engquist  and  Majda  [9], 
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dndr  dr^ 


1  \ 
2dT^j 


w  -  0 


(58) 


where  w  is  a  field  quantity  which  is  tangential  to  the 
absorbing  boundary,  n  is  the  normal  direction,  T  is 
the  tangential  direction,  and  t  is  time  normalized 
with  respect  to  the  speed  of  light.  The  second-order 
ABC  works  very  well  for  waves  which  are  at  or  near 
normal  incidence,  and  not  as  well  for  waves  which 
are  incident  at  grazing  angles.  This  second-order 
ABC  can  be  applied  to  all  the  edges  of  the  compu¬ 
tational  domain,  except  in  the  vicinity  of  any  media 
other  than  free  space  which  extends  beyond  the 
computational  domain. 

At  the  outer  boundary  and  within  or  next  to  a 
media  other  than  free  space,  a  first-order  ABC  is 
used,  since  the  fields  may  not  be  continuous  across 
the  material  interfaces.  This  program  was  intended 
for  use  with  infinite  perfectly  conducting  planes, 
where  the  tangential  electric  or  magnetic  fields  are 
discontinuous.  The  presence  of  dielectric  or  mag¬ 
netic  materials  at  the  outer  boundary  may  not  be 
adequately  treated  with  the  use  of  this  first-order 
ABC.  The  first-order  ABC  is  characterized  by  the 
following  equation. 


where  n  is  the  normal  direction.  In  order  to  simplify 
the  implementation  of  the  first-order  absorbing 
boundary.  Equation  (59)  is  differentiated  with  re¬ 
spect  to  time,  yielding 

Equation  (60)  is  equal  to  Equation  (58)  without 
the  tangential  components. 

At  the  corners  of  the  computational  domain,  the 
second-order  ABC  cannot  be  used,  since  the  normal 
and  tangential  directions  are  not  defined.  In  treating 
the  corners,  an  average  of  the  two  limiting  cases  for 
the  normal  direction  is  used  with  a  first-order  ABC, 


1 
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/_d_  /jl  +  AV 

\(9ni  dr/  \d«2  dr/ 


vv  =  0 


(61) 


the  edges,  because  there  is  only  one  point  at  each 
corner,  while  there  are  many  points  on  the  edges. 
Thus,  in  general,  using  the  first-order  ABC  at  the 
comers  will  not  cause  significant  deterioration  in 
the  accuracy  of  the  algorithm. 

For  the  //-field  polarization,  the  absorbing 
boundary  conditions  are  applied  to  the  magnetic 
field,  //,..  Applying  center  differences  to  Equation 


( 58 )  for  the  +z  boundary  k 


,  and  then 


using  temporal  and  spatial  averages  to  obtain  Hy  at 
positions  and  times  where  Hy  is  calculated,  yields 
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where  /ti ,  and  «2  are  the  two  limiting  cases  for  the 
normal  direction  at  the  corner.  This  equation  is  also 
differentiated  with  respect  to  time  so  that  the  equa¬ 
tion  has  the  same  form  as  the  second-order  ABC. 
The  ABC  at  the  comers  is  not  as  important  as  at 


The  difference  equations  satisfying  the  ABC  for  the 
other  three  edges  of  the  computational  domain  have 
the  same  form,  and  hence,  are  not  shown  here. 
The  difference  equation  satisfying  the  ABC 
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[Equation  (61)]  at  the  +a-,  +z  corner  ^ 
+  -  ,  A:  =  +  -  j  ,  is 

+  ^LZA  ^  _  I] 

At  +  a[  •  \  2  2} 


(,  =  r 


+  \,K* 

+lj 


+  //r''-(/  +  ^,A^  +  ^j 


The  absorbing  boundary  is  not  required  to  absorb 
the  incident  plane  wave,  because  the  scattered  fields, 
from  which  the  incident  plane  wave  has  been  ana¬ 
lytically  removed,  are  calculated  outside  the  inner 
region.  With  infinite  ground  planes,  the  reflected 
plane  wave  is  also  analytically  removed  from  the 
scattered  field  region  and  need  not  be  absorbed  by 
the  absorbing  boundary.  Since  the  scattered  field 
radiates  from  the  scatterer,  the  outgoing  scattered 
fields  will  not  be  incident  on  the  absorbing  boundary 
at  grazing  angles,  which  allows  the  ABC  to  absorb 
a  large  percentage  of  the  outgoing  scattered  fields. 

2.7.  Limitations 


+ 


Ar  -  A 
Ar  +  A 


Again,  the  difference  equations  at  the  other  three 
corners  have  the  same  form  and  are  not  shown. 

A  first-order  ABC  applied  at  the  +z  edge  = 

K*  4-  j  near  a  media  which  extends  beyond  the 
computational  domain  would  have  the  form 


While  the  FD-TD  technique  can  provide  very  ac¬ 
curate  predictions  and  solutions  to  various  electro¬ 
magnetic  phenomena,  there  are  a  number  of  ap¬ 
proximations  inherent  to  the  technique.  The  center 
difference  approximations  applied  to  Maxwell's 
equations  are  accurate  to  the  second  order  [i.e.,  the 
error  term  is  proportional  to  the  square  of  the  grid 
size  (A^)].  This  error  can  be  kept  to  a  minimum 
by  choosing  the  grid  size  to  be  less  than  a  tenth  of 
the  shortest  wavelength  of  interest.  Some  dispersion 
is  introduced  by  the  discretization  of  Maxwell's 
equations  (i.e..  different  frequency  components 
travel  at  slightly  different  velocities).  Geometries  are 
discretized  on  a  rectangular  grid.  Hence,  the  di¬ 
mensions  of  any  scatterer  are  restricted  to  integer 
multiples  of  the  grid  size,  and  curv  ed  surfaces  must 
be  approximated  with  staircases.  In  general,  as  long 
as  the  grid  size  is  small  compared  to  the  dimensions 
of  the  scatterer,  the  scatterer  can  be  adequately 
modeled  on  rectangular  grids.  In  order  to  obtain 
absolute  field  quantities  when  line  current  sources 
are  used,  the  sources  must  be  calibrated.  When  a 
line  source  is  too  close  to  a  boundary  ,  inaccuracies 
may  occur  in  the  field  quantities  due  to  the  inade¬ 
quate  modeling  of  the  reactive  fields.  At  the  outer 
boundary  some  reflections  will  occur  because  the 
absorbing  boundary  condition  is  an  approximate 
condition.  These  reflections  can  be  kept  tolerable 
by  locating  the  outer  boundary  at  least  half  a  wave¬ 
length  away  from  the  scatterer  [  8  ] . 
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3.  CALCULATION  OF  RADIATION/ 
SCATTERING  PATTERNS 

The  calculation  of  the  radiation  or  scattering  pat¬ 
terns  requires  the  use  of  Huygens’  principle.  Huy¬ 
gens’  principle  states  that  the  field  solution  in  a  re¬ 
gion  is  completely  determined  by  the  tangential 
fields  on  a  surface  which  completely  surrounds  the 
region  [13].  The  mathematical  formulation  of 
Huygens’  principle  for  free  space  in  two  dimensions, 
assuming  an  e~‘“'  time  dependence,  has  the  follow¬ 
ing  forms: 

£(p)=f  dS'{iuitiog(p,p')[nXH(p')] 

(65) 

+  ^  X  g(p,  P  ')[«  X  E{p  ')]} 
H(p)=  f  dS'{-/w,og(p,p')lnX  £(p')] 

(66) 

+  VXg(p,p')[nX/?(p')]} 
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where 


where: 


^(p,p')  =  -//i''(A:o|p-p'|)  (67) 

p  =  .VA  -h  fr  (68) 

p'=xx'  +  ::'  (69) 

ko  =  aiV<oPo  (70) 

and  n  is  the  normal  to  the  surface  S',  and  //o'*  is 
the  zeroth-order  Hankel  function  of  the  first  kind. 

Equations  (65)  and  (66)  are  duals,  so  that  by 
solving  one,  the  other  can  be  found  by  duality.  As¬ 
suming  that  the  integration  surface  is  a  rectangular 
box  of  dimensions  2a  in  the  x  direction  and  2b  in 
the  r direction.  Equation  (65)  becomes  the  following 
equation: 

f(p)  = -.V  j  voHA'p') 
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- //r(r)  =  -//!”(  n  (73) 
//;.."(n  =  /-,(f) +  /)',„(  f)  (74) 

The  zeroth-  and  first-order  Bessel  and  Neumann 
functions  are  approximated  by  the  following  equa¬ 
tions  [14]: 


1 


for  -3  f  <  3 


for  0  <  f  s  3 


where: 


ao(  I )  =  1 .0,  ao(  2 )  =  -2.2499997, 
ao(3)  =  1.2656208,  ao(4)  =  -0.3163866, 
ao(5)  =  0.0444479,  ao(6)  =  -0.0039444, 
tZo(7)  =  0.0002100 

6o(l)  =  0.36746691,  bo{2)  =  0.60559366, 
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bo{3)  =  -0.743200384, 

bu(4)  =  0.253001 17,  boiS)  =  -0.04261214, 

bo{6)  =  0.00427916, 
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6,(3)  =  2.1682709, 

6,(4)  =  -1.3164827,  6,(5)  =  0.3123951, 
6,(6)  =  -0.0400976,  6,(7)  =  0.0027873 
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for  3  <  f  <  X 


for  3  ^  f  <  X 

/o(  1 )  =  0.79788456,  /o(2)  =  -0.00000077, 
/n(3)  =  -0.00552740, 

^(4)  =  -0.00009512.  /o(5)  =  0.00137237, 
/o(6)  =  -0.00072805,  /o(7)  =  0.00014476 
0o(  1 )  =  -0.78539816,  00(2)  =  -0.04166397, 
0o(3)  =  -0.00003954, 

0o(4)  =  0.00262572,  0o(5)  =  -0.00054125. 
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In  the  finite  difference-time  domain  scheme,  it  is 
relatively  simple  to  obtain  the  fields  over  a  closed 
surface.  For  a  finite  geometry,  the  scattered  fields 
can  be  sampled  on  a  rectangular  box  which  encloses 
the  entire  geometry  (Figure  6).  For  an  infinite 
ground  plane  geometry',  the  scattered  fields  are  sam¬ 
pled  on  three  sides  of  a  box  which  encloses  any  dis¬ 
continuities  in  oi  above  the  infinite  ground  plane 
(Figure  7).  Image  theory  allows  the  replacement  of 
the  infinite  ground  plane  with  image  sources.  The 
appropriate  image  fields  are  shown  for  infinite  per¬ 
fect  electric  and  magnetic  conducting  planes  in  Fig¬ 
ure  9.  From  Figure  9,  it  is  apparent  that  the  tangen¬ 
tial  fields  have  been  obtained  on  a  closed  surface. 
Note  that  the  field  solution  obtained  using  the  image 
sources  is  not  valid  in  the  image  half  space. 

Sinusoidal  and  Gaussian  pulse  time-dependent 
excitations  are  treated  differently  in  order  to  obtain 
the  necessary  time  harmonic  complex  amplitudes 
of  the  fields  on  a  closed  surface.  For  sinusoidal  time 
dependent  excitations,  it  is  relatively  simple  to  ob- 
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Figure  9  Image  Sources  for  the  f-Field  Polarization  With  an  Infinite  Ground  Plane,  (a)  Electric 
Conductor,  (b)  Magnetic  Conductor. 


tain  the  complex  amplitudes  of  the  fields  by  sam¬ 
pling  the  fields  after  steady  state  has  been  reached. 
The  amplitudes  can  be  obtained  by  recording  the 
maximum  values  of  the  fields.  The  relative  phases 
can  be  obtained  by  recording  and  comparing  the 
time  of  the  maximum  values  of  the  fields.  For 
Gaussian  pulse  time-dependent  excitations,  Fourier 
transformation  of  the  excitation  source  and  the  fields 
is  performed.  At  a  particular  frequency,  by  dividing 
the  complex  Fourier  amplitudes  of  the  fields  by  the 
Fourier  amplitude  of  the  excitation,  the  complex 
amplitudes  of  the  fields  can  be  obtained  for  an  ex¬ 
citation  source  with  unity  amplitude.  For  sinusoidal 
time-dependent  excitations,  fields  and,  hence,  ra¬ 
diation  patterns  can  only  be  obtained  for  a  single 
frequency.  However,  for  Gaussian  pulse  time-de¬ 
pendent  excitations,  fields  and  radiation  patterns  can 
be  obtained  for  multiple  frequencies.  Some  limita¬ 
tions  to  the  range  of  frequencies  that  can  be  analyzed 
are  that  the  grid  must  be  fine  enough  to  adequately 
model  the  frequency  of  interest  and  that  the  Gauss¬ 
ian  pulse  must  contain  a  significant  amount  of  en¬ 
ergy  at  that  frequency. 

4.  NUMERICAL  RESULTS  AND 
DISCUSSION 

Various  electromagnetic  phenomena  were  examined 
and  animated  utilizing  the  FD-TD  code.  A  few  of 
the  electromagnetic  phenomena  which  are  exam¬ 
ined  using  this  code  include:  leading  and  trailing 
edge  diffractions  from  a  conducting  strip  and  their 
dependence  on  polarization,  creeping  waves  around 
a  conducting  cylinder,  slit  diffraction  with  plane- 
wave  excitation,  propagation  through  dielectric  and 
lossy  media,  excitation  of  a  waveguide  by  a  line 


source,  and  interaction  of  multiple  line  current 
sources.  In  each  of  the  surface  plots,  shown  in  Fig¬ 
ures  10-16,  the  outlines  of  the  scatterers  and  the 
line  sources  will  be  represented  by  a  uniform  small 
positive  height.  The  size  of  the  computational  do¬ 
main  used  to  perform  these  simulations  was  200 
nodes  X  200  nodes.  The  physical  size  of  the  com¬ 
putational  domain  was  chosen  to  be  10  meters 
X  10  meters.  Hence,  Ax  -  0.05  m  and  Az 
=  0.05  m. 

The  scatterer  in  the  first  two  cases  is  a  thin  perfect 
electric  conducting  strip  which  is  in  the  center  of 
the  computational  domain.  The  strip  is  3.5  meters 
wide  and  0. 1  meters  thick.  The  excitation  source  is 
a  Gaussian  pulse  plane  wave  which  is  incident  at  7 
degrees  above  the  plane  of  the  strip.  The  pulse  width 
is  0.354  meters  and  is  defined  to  be  equal  to  two 
standard  deviations  of  the  Gaussian  pulse.  The  shape 


Strip 

Figure  10  Leading  Edge  Diffraction  from  a  Perfectly 
Conducting  Strip  for  E-Field  Polarization. 
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of  ihc  Gaussian  pulse  is  given  by  ihc  following 
equation 

/(,v)  =  c- (84) 

where  .Vo  =  0.177  m.  The  first  case  involves  the  /;- 
field  polarization  and  is  illustrated  in  Figure  10.  For 
Figure  10,  the  electric  field,  £, ,  which  is  in  the  uni¬ 
form  and  infinite  direction,  is  plotted.  The  Gaussian 
pulse  plane  wave  appears  as  a  fairly  uniform  hill 
which  stretches  across  the  entire  computational  do¬ 
main.  The  pulse  is  moving  from  left  to  right  on  the 
figure.  The  surface  plot  was  recorded  when  the 
Gaussian  pulse  plane  wave  had  already  propagated 
near  the  end  of  the  strip.  The  ring,  which  is  centered 
at  one  edge  of  the  strip,  is  the  leading  edge  diffraction 
from  the  edge  of  the  strip  which  was  illuminated 
first  by  the  plane  wave.  The  diffraction  from  the 
leading  edge,  in  this  case,  is  strong  because  there  is 
a  strong  discontinuity  in  the  electric  field  w'hcn  the 
pulse  is  incident  on  the  leading  edge.  The  perfect 
electric  conductor  requires  the  tangential  electric 
field  at  the  surface  to  be  zero.  Hence,  when  the  in¬ 
cident  electric  field,  which  is  parallel  to  the  strip 
surface,  impinges  on  the  strip,  a  strong  reflected  wave 
is  generated  in  order  to  satisfy  the  boundary  con¬ 
ditions.  At  the  trailing  edge,  the  boundary  condition 
states  that  the  current  must  be  continuous.  The  sur¬ 
face  current.  A',  is  defined  to  be 

K=nXH  (85) 

Since  the  magnetic  field  is  almost  perpendicular  to 
the  strip,  the  tangential  magnetic  fields  at  the  surface 
will  be  weak.  Thus,  the  surface  currents  will  also  be 


Strip  ^ 

Figure  II  Trailing  Edge  Diffraction  From  a  Perfectly 
Conducting  Strip  for  //-Field  Polarization. 


F'igure  12  Creeping  Waves  Around  a  Perfectly  Con¬ 
ducting  Cylinder  for  //-Field  Polarization. 


weak,  and  there  will  be  little  diffraction  at  the  trailing 
edge.  The  second  case,  which  is  illustrated  in  Figure 
1 1.  involves  the  //-field  polarization.  In  this  case, 
the  magnetic  field,  //, .  is  plotted.  The  Gaussian  pulse 
plane  wave  has  propagated  past  the  strip,  and  now 
two  circular  rings  can  be  seen.  The  larger  ring  is  the 
leading  edge  diffraction,  which  is  weaker  than  the 
prex  ious  case  because  the  electric  field  of  the  incident 
plane  wave  points  nearly  perpendicularly  to  the 
strip,  and  causes  a  less  severe  discontinuity  in  the 
tangential  electric  field  at  the  strip  surface.  The 
smaller  ring  represents  the  trailing  edge  diffraction. 
In  this  case,  the  surface  currents  will  be  particularly 
strong,  because  the  magnetic  field  is  always  parallel 
to  the  strip’s  surface.  The  strong  trailing  edge  dif¬ 
fraction  is  due  to  the  large  discontinuity  in  the  sur¬ 
face  current  at  the  trailing  edge  of  the  strip. 

The  third  scatterer  involves  a  perfectly  conduct¬ 
ing  cylinder.  The  excitation  source  is  again  a  Gauss¬ 
ian  pulse  plane  wa\e  which  is  incident  at  90  degrees 
for  the  //-field  polarization.  The  cylinder  has  a  di¬ 
ameter  of  2.5  meters,  and  the  pulse  width  is  again 
0.354  meters.  Figure  1 2  illustrates  the  magnetic  field. 
//, ,  as  the  plane  wave  scatters  off  the  cylinder.  Again, 
the  incident  plane  wave  can  be  seen  after  it  has 
propagated  past  the  cylinder.  The  circular  arc  is 
traveling  outward  away  from  the  cylinder  and  is  the 
specular  reflection  from  the  cylinder.  In  the  vicinity 
of  the  shadow  region,  where  the  cylinder  blocks  the 
incident  plane  w-ave,  the  circular  arc  is  still  attached 
to  the  conducting  cylinder  and  travels  around  the 
back  of  the  cylinder.  This  phenomenon  corresponds 
to  a  creeping  wave.  Creeping  waves  occur  much 
more  strongly  for  the  //-field  polarization,  because 
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Figure  13  Diffraction  of  a  Plane  Wave  Through  a  Slit 
in  an  Infinite  Perfectly  Conducting  Ground  Plane  for  //- 
Field  Polarization. 


Stronger  surface  currents  are  generated,  which  es¬ 
sentially  allow  the  fields  to  remain  attached  to  the 
cylinder. 

The  fourth  case  illustrates  diffraction  of  a  plane 
wave  through  a  slit  in  an  infinite  perfect  electric 
conducting  ground  plane.  The  infinite  ground  plane 
is  O.l  meters  thick  and  the  slit  is  0.3  meters  wide. 
The  plane  wave  is  normally  incident  upon  the 
ground  plane  from  the  upper  half  space  and  has  a 
Gaussian  pulse  time  dependence  with  a  pulse  width 
of  0.4  meters.  In  this  case,  the  //-field  polarization 
is  examined,  and  Figure  1 3  shows  the  magnetic  field, 
//,..  The  reflected  plane  wave  is  evident  and  has  only 
been  slightly  perturbed  by  the  presence  of  the  slit. 
The  fields  that  pass  through  the  slit  radiate  cylin- 
drically  from  the  aperture,  thus  illustrating  slit  dif¬ 
fraction. 

Figure  14  illustrates  the  propagation  of  a  Gauss¬ 
ian  pulse  plane  wave  through  a  dielectric  cylinder 
for  the  //-field  polarization.  The  cylinder  has  a  di¬ 
ameter  of  2.5  meters,  and  the  pulse  width  is  again 
0.354  meters.  The  dielectric  cylinder  has  (r  =  2,  n, 
=  1 ,  and  a  =  a„,  =  0.  Figure  14  shows  the  magnetic 
field,  Hy.  The  specular  reflection  off  the  cylinder  is 
not  as  strong  as  the  specular  reflection  for  the  per¬ 
fectly  conducting  cylinder  because  some  of  the  fields 
are  reflected  and  some  are  transmitted.  The  fields 
within  the  dielectric  cylinder  lag  behind  the  fields 
in  free  space  because  they  propagate  at  a  slower  ve¬ 
locity  in  the  dielectric  than  in  free  space.  The  mag¬ 
netic  fields  within  the  cylinder  have  higher  intensities 
for  the  following  reasons.  First,  the  phase  matching 
condition  on  the  surface  of  the  dielectric  cylinder 


tends  to  bend  the  electromagnetic  waves  toward  the 
bisector  of  the  cylinder  which  is  parallel  to  the  in¬ 
cident  wave  vector,  which  focuses  the  fields  to  some 
extent.  The  reflection  coefficient  at  a  planar  bound¬ 
ary  between  free  space  and  a  dielectric  with  relative 
permittivity  greater  than  1  for  the  //-field  polariza¬ 
tion  is  always  greater  than  zero.  Hence,  the  trans¬ 
mission  coefficient  will  always  be  greater  than  1, 
and  the  magnetic  fields  transmitted  through  the  di¬ 
electric  will  be  greater  than  the  incident  magnetic 
fields.  Treating  each  point  of  the  dielectric  cylinder 
locally  as  a  planar  interface,  it  is  clear  that  the  mag¬ 
netic  fields  within  the  cylinder  will  be  greater  than 
the  incident  fields  outside  the  cylinder. 

The  next  case  examines  the  excitation  of  the  per¬ 
fectly  conducting  parallel  plate  waveguide  by  a  line 
current  source.  The  waveguide  is  0.55  meters  tall 
and  assumed  to  be  infinite  in  length,  and  the  thick¬ 
ness  of  the  walls  is  0. 1  meters.  The  line  source  is 
vertically  centered  within  the  cavity  and  is  located 
in  the  middle  of  the  section  of  the  waveguide  shown. 
The  //-field  polarization  is  examined  and  magnetic 
field.  Hy,  is  plotted  in  Figure  15.  The  line  current 
source  radiates  at  a  frequency  of  200  MHz.  The 
wavelength  in  free  space  is  1.5  meters.  The  line  cur¬ 
rent  source  excites  the  TEM  mode  within  the  wave¬ 
guide.  Two  waves,  which  travel  in  opposite  direc¬ 
tions,  are  launched  by  the  line  source.  There  is  also 
a  lossy  dielectric  material  which  spans  the  height  of 
the  waveguide  and  is  0.9  meters  wide.  The  lossy 
material  has  i,-tir=  1,  ff  =  0.1,  and  a„  =  0.  The 
fields  are  attenuated  as  they  propagate  through  the 
dielectric  and  it  is  apparent  that  the  fields  that  have 
passed  through  lossy  dielectric  are  smaller  than  the 
fields  which  have  not. 

The  last  case  illustrates  the  interference  pattern 
generated  by  four  line  sources.  Figure  16  plots  the 


Figure  14  Propagation  Through  a  Dielectric  Cylinder 
With  (,=  2,  fjr  =  I,  and  =  am  =  0  for  //-Field  Polar¬ 
ization. 
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Figure  IS  Excitation  of  a  Parallel  Plate  Waveguide  by 
a  Line  Current  Source,  and  Propagation  Through  a  Lossy 
Dielectric  With  t,  =  n,  =  1 ,  o-  =  0. 1 ,  and  a,„  =  0  for  H- 
Field  Polarization. 


electric  field,  Ey,  for  the  £-field  polarization.  The 
line  sources  are  excited  at  a  frequency  of 333. 3  MHz. 
In  free  space,  the  wavelength  is  0.9  meters.  The  line 
current  sources  are  separated  by  0.45  meters,  or  half 
a  wavelength.  The  line  sources  are  in  phase  and  have 
equal  magnitudes  or.  in  other  words,  have  the  exact 
same  time  dependence.  On  the  line  on  which  the 
line  sources  lie,  the  radiation  is  very  small  since  the 
radiation  from  each  line  source  almost  exactly  can¬ 
cels  the  radiation  from  the  adjacent  line  sources.  On 
this  line,  the  radiation  from  each  line  source  will  be 
1 80'  out  of  phase  with  the  next  line  source,  because 
they  are  each  spatially  separated  by  half  a  wavelength 
in  this  direction.  On  the  line  perpendicular  to  the 
line  connecting  the  four  line  sources,  the  radiation 
from  the  line  sources  will  add  constructively.  The 
waves  reinforce  each  other  because  they  are  in  phase 
since  every  point  on  this  line  is  almost  equidistant 
to  the  four  line  sources.  Hence,  the  radiation  will 


Figure  16  Interference  Pattern  Produced  by  Four  Line 
Sources  for  £-Field  Polarization. 


be  largest  in  these  two  directions.  There  will  also  be 
additional  nulls  and  maxima  as  the  four  line  sources 
add  constructively  and  destructively  for  different 
directions. 

Due  to  the  limitations  of  the  surface  plots  and 
the  inability  to  include  color  plots,  it  was  not  possible 
to  show  all  the  capabilities  of  the  code  in  this  study. 
One  extremely  useful  aspect  of  this  code  is  the  ca¬ 
pability  to  calculate  and  successively  display  the 
fields  in  color  in  order  to  visually  produce  the  prop¬ 
agation  of  electromagnetic  waves  in  real  time.  The 
typical  computational  time  for  animation  consisting 
of  500  time  steps  and  100  color  plots  is  less  than  1 
hour  on  a  33-MHz  IBM  386  compatible  personal 
computer.  After  the  FD-TD  compulation  has  been 
completed,  the  actual  animation  can  be  played  back 
at  about  five  frames  per  second. 

In  addition  to  the  capability  to  visualize  the  fields, 
the  computer  code  can  generate  radiation  and  scat¬ 
tering  patterns.  Figures  17-20  illustrate  radiation 
from  a  linear  array,  and  scattering  from  a  perfectly 
conducting  strip  and  a  large  slit  in  an  infinite  ground 
plane.  In  each  of  these  polar  plots,  the  field  ampli¬ 
tude  is  represented  in  dB  for  a  radius  of  50  meters. 
For  £-field  polarization  cases,  £,  is  plotted,  and  a 
complex  amplitude  of  1  V/m  is  assumed  for  the 
excitation  source.  For  //-field  polarization  cases,  //, 
is  plotted,  and  a  complex  amplitude  of  1  A/m  is 
assumed  for  the  excitation  source. 

The  radiation  pattern  of  the  four-element  array 
previously  examined  in  Figure  16  is  illustrated  in 
Figure  1 7.  Four  line  current  sources  for  the  £-field 
polarization  are  excited  at  a  frequency  of  333.3 


Figure  17  Radiation  Pattern  of  a  Four-Element  Array 
of  Line  Sources  for  £-Field  Polarization 
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MHz.  and  the  sources  are  separated  by  0.45  meters 
or  half  a  wavelength.  The  line  sources  have  the  same 
phase  and  magnitude.  The  radiation  pattern  has 
four-fold  symmetry  since  the  problem  has  four-fold 
symmetry.  The  radiation  pattern  has  nulls  in  the 
+A-,  -,v  directions,  60“  and  120“  above  and  below 
the  jc-axis.  The  main  lobes  are  in  the  -hr  and  -s 
directions.  The  main  lobes  correspond  to  a  electric 
field  of 0.3096  V/m.  Recall  that  the  line  sources  are 
uncalibrated,  so  that  the  absolute  numbers  will  be 
dependent  on  the  grid  size.  The  radiation  pattern 
obtained  here  is  consistent  with  the  radiation  pattern 
that  would  be  obtained  using  linear  array  theory. 

Figure  1 8  shows  the  scattering  pattern  of  perfect 
electric  conducting  strip.  The  strip  is  illuminated  by 
a  300-MHz  sinusoidal  plane  wave  at  an  angle  of  45“ 
above  the  x-axis  for  the  £-field  polarization.  The 
dimensions  of  the  strip  are  3.1  m  X  0.1  m.  As  ex¬ 
pected,  the  main  lobes  of  the  scattering  pattern  occur 
in  the  specular  direction  ( 135“  above  the  .v-axis) 
and  the  forward  scatter  direction  ( 135“  below  the 
.v-axis).  The  maximum  scattered  electric  field  is 
0.268  V/m. 

Figures  19  and  20  show  the  radiation  patterns 
for  a  large  slit  illuminated  by  a  Gaussian  pulse  plane 
wave  at  normal  incidence  from  the  lower  half  space 
for  the  //-field  polarization.  The  slit  is  4.9  m  wide 
and  the  infinite  perfect  electric  conducting  plane  is 
0.1  m  thick.  Since  a  Gaussian  pulse  time  dependence 
is  used,  multiple  frequencies  can  be  analyzed.  The 
frequencies  examined  in  Figures  19  and  20  are  300 
MHz  and  400  MHz,  respectively,  and  the  corre- 
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Figure  18  Scattering  Pattern  of  a  Perfectly  Conducting 
Strip  Illuminated  at  a  45“  Grazing  Angle  for  £-Field  Po¬ 
larization. 


Figure  19  Diffraction  Pattern  of  a  Slit  in  an  Infinite 
Ground  Plane  Illuminated  at  Normal  Incidence  at  300 
MHz  for  //-Field  Polarization. 


spending  wavelengths  are  1  m  and  0.75  m.  The 
maximum  peaks  in  the  radiation  pattern  occur  in 
+z  direction  and  are  0.650  A/m  at  300  MHz  and 
0.809  A/m  at  400  MHz.  Fraunhofer  diffraction 
states  that  the  fields  in  an  aperture  are  related  to  the 
far-fieJd  pattern  by  a  Fourier  transform  [  1 3 ) .  In  these 
cases,  the  field  distribution  is  essentially  uniform 
over  the  aperture,  so  the  far-field  patterns  should 


Figure  20  Diffraction  Pattern  of  a  Slit  in  an  Infinite 
Ground  Plane  Illuminated  at  Normal  Incidence  at  400 
MHz  for  //-Field  Polarization. 
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the  aperture  is  electrically  larger,  so  its  transform, 
and  hence,  radiation  pattern,  is  a  narrower  sine 
function. 

5.  SUMMARY 

This  FD-TD  code  can  be  used  to  examine  various 
electromagnetic  phenomena.  The  code  successfully 
and  accurately  predicted  the  interaction  between 
various  sources  and  objects.  Because  the  FD-TD 
technique  has  the  capability  to  animate  propagation, 
diffraction  and  scattering  of  electromagnetic  waves, 
it  is  well  suited  as  an  educational  tool  for  students 
of  electromagnetic  wave  theory  and  its  applications. 
The  additional  capability  of  generating  radiation 
patterns  provides  an  alternative,  more  quantitative 
way  of  viewing  electromagnetic  radiation. 
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Time  Domain  Modeling  of  Impedance  Boundary 
Condition 


oiquc  is  discussed  and  the  overall  scheme  is  verified  numerically 
for  a  one-diniensional  configuration. 
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Abstract— A  methodology  developed  to  handle  dispersive  materials  in 
the  time  domain  is  c.xtended  to  model  the  dispersive  characteristics  of 
the  impedance  boundary  condition  used  for  a  thin  layer  coating  over 
perfect  conductors.  The  impedance  boundary  condition  is  first  approx¬ 
imated  as  a  rational  function  of  frequency.  This  rational  function  is 
then  transformed  to  a  time  domain  equation,  resulting  in  a  partial 
differential  equation  in  space  and  time.  Discretization  of  the  time  do¬ 
main  model  to  efficiently  handle  the  thin  layer  coating  is  presented  in 
the  context  of  the  finite-difference  time-domain  (FD-TD)  technique.  The 
methodology  is  verified  by  solving  a  one-dimensional  problem  using  the 
FD-TD  technique  and  comparing  with  the  analytical  results. 


I.  Introduction 

Electrically  fine  structures  often  appear  in  practical  applications. 
To  resolve  the  electromagnetic  behavior  of  these  structures,  very 
fine  grids  are  needed  in  numerical  techniques  (e.g..  the  finite-dif¬ 
ference  time-domain  technique).  Alternatively,  one  may  incorpo¬ 
rate  the  localized  physical  behavior,  such  as  the  impedance  bound¬ 
ary  condition  [1]  and  thin  wire  formulations  [2],  [3],  of  these  fine 
structures  into  discretization  schemes. 

Thin  surface  coatings  on  metallic  bodies  appear  in  many  .scatter¬ 
ing  problems.  In  principle,  these  thin  surface  coatings  can  be  mod¬ 
eled  numerically  and  geometrically  by  very  fine  grids  with  appro¬ 
priate  discretization  schemes.  The  disadvantage  associated  with 
such  an  approach  is  the  large  computer  memory  requirement.  Fur¬ 
thermore.  in  the  finite-difference  time-domain  (FD-TD)  technique, 
the  time  increment  is  usually  determined  by  the  smallest  grid  size 
in  the  entire  computational  domain  to  satisfy  the  stability  condi¬ 
tion. 

In  modeling  a  thin  layer  coating  in  the  frequency  domain,  the 
concept  of  impedance  boundary  condition  can  be  used  to  avoid  the 
fine  layers  of  grids.  The  impedance  boundary  condition  relates  the 
tangential  fields  on  the  coating  to  their  normal  derivatives,  which 
is  derived  using  the  configuration  of  a  half-space  conductor  with 
thin  layer  coating.  The  resulting  impedance  boundary  condition  is 
frequency  dependent,  and  this  dispersive  nature  of  the  impedance 
boundary'  condition  causes  difficulty  in  the  time  domain  modeling. 
In  this  paper,  a  time  domain  technique  used  to  treat  dispersive  ma¬ 
terials  is  employed  to  convert  the  impedance  boundary  condition 
to  the  time  domain  |4].  Following  the  idea  in  (5|.  the  impedance 
boundary  condition  is  approximated  by  a  rational  function  of  the 
frequency.  Then,  the  tangential  fields  arc  related  to  their  normal 
derivatives  by  a  partial  differential  equation  in  space  and  time.  A 
numerical  discretization  scheme  in  the  context  of  the  FD-TD  tcch- 
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II.  Time  Domain  Description  oe  the  I.mpedance  Boundary 
Condition 

The  impedance  boundary  condition  of  the  coated  conductor  can 
be  derived  based  on  a  two-layer  configuration  (Fig.  1).  By  ignoring 
the  variation  along  the  tangential  directions,  the  tangential  electric 
field  can  be  related  to  its  no.mal  derivative  by  the  following  equa¬ 
tion: 


tir  ! - 

E  =  —  tan  «:  As/m^) ,  (1) 

k  oy 

where  k  is  the  free  space  wavenumber,  n,  and  e,  are  the  relative 
permeability  and  relative  permittivity,  respectively,  and  ij,  = 
'In,! t,  is  the  relative  impedance  (impedance  normalized  to  rjo). 
Inverse  Fourier  transformation  may  be  used  to  convert  the  above 
equation  to  the  lime  domain.  However,  it  is  relatively  complicated 
and  the  result  may  not  be  suitable  for  numerical  analysis.  Instead, 
following  the  procedure  outlined  in  [Sj.  the  above  equation  is  ap¬ 
proximated  using  a  rational  function  of  the  frequency.  With  the 
substitution  of  -to.'  by  d/d/,  the  time  domain  description  of  the 
impedance  boundary  condition  is  obtained. 

The  rational  function  approximated  of  (1)  can  be  obtained  by 
expressing  the  tangent  function  as  the  ratio  of  sine  and  cosine.  Next, 
the  Taylor  series  expansions  of  these  two  functions  are  used  to 
obtain  the' rational  function  form.  By  keeping  the  first  two  terms  in 
the  Tay  lor  series  expansions  of  the  sine  and  cosine  functions,  the 
first-order  rational  function  approximation  is  obtained. 


E  =  »),  A  sVT, 


'  ~6^'-^'^'‘-dE 
1  -  ];k'  A-m,«, 


(2) 


Similarly,  by  keeping  the  first  three  terms  of  Taylor  series  expan¬ 
sion.  the  second-order  approximation  can  be  obtained: 


E  =  9,  A 


k  A 


1  ’  /  -’  X  J  4.  '  1  ^  \  ^ 


' 

T--  (3) 


In  general,  both  the  relative  pemiittivity  and  permeability  in  (2) 
and  (3)  can  be  complex  to  account  for  electric  and  magnetic  losses. 
However,  in  this  paper,  we  will  consider  only  the  loss  due  to  elec¬ 
trical  conductivity. 

The  time  domain  expressions  corresponding  to  (2)  and  (3)  can 
be  obtained  by  substituting  —iie  by  d  idi.  or  equivalently,  —ik  by 
d/dr.  where  t  =  r,)/  is  the  normalized  time.  .Assuming  electrical 
conduction  loss  only,  we  obtain 


Ir 


.  otj„ 

+  'T- 


(4) 


where  t',  a  and  ijo  ^re  dielectric  constant,  electric  conductivity  and 
free  space  impedance,  respectively.  The  first-order  time  domain 
description  of  the  impedance  boundary  condition  corresponding  to 
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Fig.  1.  Two-layer  configuration. 


(2)  is  given  by 


1  ,31,  3- 

1  + 


=  A 


.  >  .2  3  1.2,  3- 

1  + 


3v 


(5) 


Applying  the  same  transformation  to  (3),  the  second-order  time 
domain  description  of  the  impedance  boundary  condition  is  given 
by 


1.4.  o  I  .  ■>  ,1  d' 


3r- 


=  A 


1  -R 


1  ,  3  /I  ,  I  ,  ,  ,N  3- 

gA-o,o^  4-(^gA-e  +n5Ao-,5j^ 


'  4  ,  5  1  .  4  3" 

+  -  €  OTin  +  -  A  <  “ 

60  ’“St'  120  dr* 


3£ 

3v 


(6) 


The  above  time  domain  descriptions  of  the  impedance  boundary 
condition  arc  panial  differential  equations  in  space  and  time.  In 
fact,  these  equations  may  be  derived  directly  in  the  time  domain 
with  a  finite  difference  approximation  and  a  simple  averaging 
scheme  [4]  (see  Appendix).  The  impedance  boundary  condition  ac¬ 
counts  for  the  interaction  between  air-dielectric  interface  and  con¬ 
ducting  surface.  In  the  time  domain  there  is  a  time  delay  associated 
with  this  interaction.  This  delay  is  partially  modeled  by  (5)  and 
(6). 


in.  Discretization 

There  arc  many  discretization  schemes  one  can  use  to  discretize 
(5)  and  (6)  for  the  finite-difference  time-domain  technique.  The 
discretization  of  the  normal  derivative  is  quite  important  in  the  im¬ 
plementation  of  the  impedance  boundary  condition.  In  order  to  have 
an  accurate  approximation,  the  vanishing  tangential  electric  field 
on  the  conducting  surface  is  used.  This  leads  to  interpolation  of  the 
electric  field  in  three  or  more  locations.  Employing  the  Lagrange 
interpolation  formula  for  three  electric  fields  (Fig.  2),  the  normal 
derivative  is  discretized  as  follows: 


^.v 


a 

A,(l  -I-  a) 


O) 


Eo  =  0  El  E  2 

Fig.  2.  Discretization  nodes  for  normal  derivative. 


In  the  above  equation  a  represents  the  ratio  of  the  layer  thickness 
to  the  grid  size.  A /A,.  It  turns  out  that  a  better  approximation  of 
the  normal  derivative  can  be  obtained  by  including  the  effect  of  the 
dielectric  constant  of  the  coating.  Because  the  wave  velocity  in  the 
coating  is  slower  than  that  in  the  free  space  where  discretization 
applied,  the  effective  thickness  of  the  coating  should  be  AyfP. 
Therefore,  the  following  equation  is  used; 

A  Je' 

a  =  - . 

The  temporal  derivatives  in  (5)  and  (6)  may  also  be  discretized  in 
many  forms.  A  simple  second-order  center-differencing  is  used  in 
this  paper. 


IV.  Numerical  Results 

Equations  (5)  and  (6)  describe  the  time  domain  modeling  of  the 
impedance  boundary  condition  given  by  (1).  Equation  (7)  and  cen¬ 
ter-differencing  in  time  provide  a  possible  discretization  scheme. 
To  validate  the  approach  outlined  in  this  paper,  a  one-dimensional 
reflection  problem  is  simulated  using  these  models  together  with 
the  finite-difference  time-domain  method  (3).  The  reflected  wave 
is  calculated  in  the  time  domain  and  Fourier  transformed  to  obtain 
the  reflection  coefficient  as  a  function  of  frequency.  The  results  are 
compared  with  the  exact  solution  obtained  directly  in  the  frequency 
domain. 

Fig.  1  shows  the  configuration  of  the  problem.  The  layer  has 
thickness  of  0.2  cm.  The  coating  material  is  assumed  to  have  di¬ 
electric  constant  of  5  and  conductivity  of  0.2  mho/meter.  The 
computation  domain  contains  500  nodes  with  grid  size  being  0.2 
cm.  The  normalized  time  increment  A,  is  0.2  cm  to  satisfy  the 
stability  criterion  and  to  minimize  numerical  dispersion.  At  the  first 
node  of  the  computational  domain,  an  incident  Gaussian  pulse 
modulated  by  a  carrier  at  8  GHz  is  imposed.  This  Gaussian  pulse 
has  a  half-power  pulse  width  of  0.0625  nano-second,  which  cor¬ 
responds  to  a  bandwidth  of  6  GHz.  The  last  node  of  the  computa¬ 
tional  domain  is  placed  at  the  free  space/dielectric  interface  where 
the  first-order  or  the  second-order  approximation  to  the  impedance 
boundary  condition  (IBC)  is  applied  to  simulate  the  thin  layer  coat¬ 
ing. 

The  exact  solution  for  the  reflection  coefficient  is  given  by 
if),  tan  {k  A  Vp,*,)  -I-  1 

R{k)  =  - — - .  (8) 

(Tj,  tan  {k  A  v/irt,)  —  1 

where  jj,  and  f,  is  complex.  For  the  assumed  parameters,  the  mag¬ 
nitude  and  the  phase  of  this  reflection  coefficient  are  plotted  in  Fig. 
3(a)  and  (b)  (solid  curves),  respectively,  from  2  GHz  to  16  GHz. 
The  results  obtained  by  using  the  FD-TD  simulation  are  also  show  n 
in  Fig.  3(a)  and  (b).  The  “O"  and  “-i-”  curves  represent  results 
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FREQUENCY  (GHz) 

(b) 

Fig.  3.  (a)  Magnitude  of  reflection  coefficient  versus  frequency,  (b)  Phase 
of  reflection  coefficient  versus  frequency. 


0.10  O.U  0.18  0.22  0.26 

COATING  THICKNESS  (cm) 

Fig.  4.  Percentage  error  for  second-order  IBC  versus  coating  thickness. 

V.  Summary 

In  this  paper,  a  time  domain  modeling  of  the  impedance  bound¬ 
ary  condition  is  derived  and  expressed  in  terms  of  partial  differ¬ 
ential  equations  in  space  and  time.  A  possible  discretization  scheme 
which  incorporates  the  effective  thickness  of  the  layer  is  presented. 
Numerical  results  indicate  the  validity  of  the  modeling  as  well  as 
the  suggested  discretization  scheme.  Although  this  model  is  only 
verified  for  a  one-dimensional  problem,  the  generalization  to  higher 
spatial  dimensions  is  possible  since  the  equation  is  applied  only  at 
the  interface.  However,  tangential  variations  may  need  to  be  con¬ 
sidered  in  those  cases.  Furthermore,  the  concept  of  using  the  ra¬ 
tional  functions  to  approximate  the  frequency  domain  response  and 
converting  to  the  time  domain  may  be  used  to  characterize  other 
fine  structures. 

Appendix 

For  simplicity,  the  relative  permittivity  and  permeability  are  as¬ 
sumed  real.  For  a  plane  wave  normally  incident  with  waveform  g. 


obtained  using  (5)  and  (6),  respectively.  Both  the  first-  and  the 
second-order  approximations  yield  good  agreement  with  the  exact 
solution  in  the  phase  of  the  reflection  coefficient.  However,  the 
advantage  of  using  the  second-order  approximation.  (6).  is  clearly 
shown  in  the  magnitude  of  the  reflection  coefficient.  The  results 
obtained  using  the  first-order  approximation  match  the  exact  solu¬ 
tion  at  low  frequencies  with  increasing  discrepancy  at  higher  fre¬ 
quencies. 

Fig.  4  shows  the  percent  error  of  the  reflection  coefficient  mag¬ 
nitude  as  a  function  of  the  coating  thickness.  These  errors  are  for 
the  FD-TD  results  with  the  second-order  approximation  of  the  IBC 
using  the  same  grid  size  and  the  same  material  parameters  while 
varying  the  coating  thickness.  The  three  curves  represent  errors  at 
three  different  frequencies.  The  errors  for  the  half-power  frequen¬ 
cies,  5  GHz  and  1 1  GHz,  are  shown  by  ••  -p  "  and  "  x ,"  respec¬ 
tively.  The  errors  at  the  carrier  frequency  arc  shown  by  “O.” 
These  errors  are  all  well  within  1.0  percent,  and  they  increase  with 
increasing  thickness,  as  expected.  It  should  be  noted  that  the  errors 
are  relatively  small  even  when  the  thickness  of  the  coating  is  up  to 
1  /5  of  the  wavelength  inside  the  layer  as  in  the  case  of  1 1  GHz  at 
the  thickness  of  0.25  cm. 


=  «(t  +  y).  =  -g(T  +  >•).  (Al) 

the  reflected  and  transmitted  waves  are 

Er  =/(t  -  .V).  TJoffr  =/(t  -  .v);  (A2) 

E,  =  p{T  -I-  y  H-  A) 

-  p(T  -  •Jurfr  y  -  A);  (A3a) 

TJoW,  =  -  —  p{T  -p  'Ip^ry  +  '/aT^A) 

Vr 


(A3b) 


In  (A2)  and  (A3),  the  subscript^  r  and  s  denote  the  reflected  wave 
and  the  wave  inside  the  coating,  respectively.  Note  that  the  bound¬ 
ary  condition  on  the  surface  of  the  conductor  has  been  satisfied. 
Imposing  the  boundary  conditions  of  continuous  tangential  electric 
and  magnetic  fields  at  the  interface, 

E  =  g(T)  -I-  /(t)  =  piT  +  A)  -  P(t  -  'fpjr  A), 


(A4a) 
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VaH  =  s(^)  +  /(^)  = 


1 

- p(T  + 

V, 


A) 


Comments  on  “Criteria  for  the  Onset  of  Oscillation 
In  Microwave  Circuits” 


(A4b) 


In  addition,  the  total  tangential  electric  and  magnetic  fields  are  re¬ 
lated  by; 

E  =  2g(T)  -F  T/oW.  (A5) 

Applying  the  center-differencing  and  averaging  approximations: 
P(t  +  yfpTr  A)  -  p(t  -  yfii^r  A) 

~  P'{t)  +  \A^(p,t,)^^-p'"(T)  (A6a) 

p(T  -F  A)  +  p(t  -  yfp^,  A) 

~  2p{T)  +  A>,e,p''(T),  {A6b) 


where  denotes  derivative  with  respect  to  the  argument.  Sub¬ 
stituting  these  back  into  (A4), 


2g(T)  =  —  p{T)  +  2A'Jp,(.,p'(r) 
V, 


A- 


a’ 


+  -  H,i,p"{T)  +  —  (t) 


(A7) 


rtoH  = 


(2/)(t)  -F  A-p,(,p''iT)]. 
Vr 


(A8) 


Substituting  (A7)  and  (A8)  to  (A5),  £  is  related  to  the  first  and 
third  derivatives  of  p(T).  Again,  using  (A8)  and  ignoring  the  terms 
involving  derivatives  higher  than  third  order  by  assuming  slowly 
varying  field  in  the  time  scale  of  A,: 


„  .  ,  A^pUr  d'rtoH 

£=  -A,.—  +  ^ - 


Using  Maxwell's  equations; 


IT  _  .  aV;<,  d'E 

^  A/4.. 

ov  3  or  dv 


(A9) 


This  is  what  one  expects  if  one  expands,  in  Taylor  series,  the  taii- 
geni  function  in  (1)  and  converts  the  expanded  equation  to  the  time 
domain  using  the  procedure  described  earlier. 
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Robert  W.  Jackson 

The  paper  listed  above'  notes  that  the  device  reflection  coefli- 
cient,  Pa'fj).  in  the  expression, 

I/-  =  - - ;; 

I  -  ra(.s)r,(i) 

represents  the  port  reflection  coefficient  of  a  device  which  may  re¬ 
sult  in  an  unstable  circuit  only  after  connecting  it  with  a  resonator 
having  a  reflection  coefficient.  r<.(s).  This  is  an  important  condi¬ 
tion  and  is  somewhat  vague  as  worded.  In  order  to  use  the  Nyquist 
criterion  to  determine  the  stability  of  the  device-circuit  combina¬ 
tion,  Fa  must  have  no  right  half  plane  poles.  This  amounts  to  in¬ 
suring  that  the  device  does  not  oscillate  into  the  reference  imped¬ 
ance  (50  ohms  for  example).  If  F^  has  been  determined  from 
measurements,  presumably  the  device  is  not  oscillating  during  the 
measurement  and  therefore  there  are  no  right  half  plane  poles. 

In  CAD  simulations  of  possibly  unstable  circuits,  the  location  of 
Fj  poles  is  not  always  so  clear.  For  a  simple  amplifier  circuit  such 
as  the  one  described  in  the  above  referenced  paper,  one  can  assume 
no  right  half  plane  poles  in  F^  if  the  S,,  and/or  S..  coefficients  of 
the  FET  have  magnitudes  less  than  one.  To  see  this,  consider  the 
partial  circuit  formed  by  a  50  ohm  lemiination  on  port  2  and  any 
passive  termination  on  port  1.  If  |S|||  <  1,  the  input  termination 
sees  a  passive  impedance  and  therefore  the  panial  circuit  is  stable. 
Since  the  partial  circuit  is  stable,  Fj  (50  ohm  reference)  seen  look¬ 
ing  in  at  port  2  has  no  poles  in  the  right  half  plane.  If.  as  in  the 
amplifier  example' .  Fj  has  a  magnitude  greater  than  1 .  the  Nyquist 
criterion  as  described  can  then  be  applied  to  study  the  stability  ef¬ 
fects  of  various  port  2  terminations.  In  simulations  using  devices 
with  extra  feedback,  oscillators  for  example,  often  the  magnitudes 
of  S|i  and  Svi  are  both  greater  than  one  and  this  approach  breaks 
down. 

A  more  generally  applicable  use  of  the  Nyquist  stability  criterion 
has  been  known  for  years,  but  the  current  widespread  use  of  mi¬ 
crowave  CAD  makes  it  must  easier  to  apply.  As  discussed  in  the 
literature  (1],  12]  the  admittance  between  any  two  nodes  in  an  ac¬ 
tive  circuit  cannot  have  right  half  plane  zeros  if  the  circuit  is  to  be 
stable.  If  one  were  to  apply  the  Nyquist  test  to  such  an  admittance, 
the  resulting  Nyquist  locus  of  points  cannot  encircle  zero  in  a 
clockwise  sense  if  the  circuit  is  stable.  It  is  trivial  for  modem  mi¬ 
crowave  CAD  programs  to  calculate  the  necessary  admittances  vs 
frequency.  Polar  plotting  of  admittances  is  not  always  available  but 
a  quick  sketch  is  easy  to  do.  It  should  be  noted  that  the  number  of 
Nyquist  encirclements  only  gives  the  difference  between  the  num¬ 
ber  of  right  half  plane  zeros  and  poles  in  the  admittance  function. 
If,  for  example,  the  admittance  at  a  particular  node  pair  has  an 
equal  number  of  right  half  plane  poles  and  zeros,  the  Nyquist  plot 
would  not  encircle  the  origin  even  though  the  circuit  is  unstable. 
Thus  a  clockwise  encirclement  insures  instability,  but  no  encircle¬ 
ment  docs  not  insure  stability.  Since  admittances  at  various  node 
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To  understand  the  physical  meaning  of  rational  reflection  coefficients  in  inverse-scattering  theory  for  optical 
waveguide  design  [J.  Opt.  Soc.  Am.  A  6,  1206  (1989)],  we  studied  the  relationship  between  the  poles  of  the 
transverse  reflection  coefficient  and  the  modes  in  inhomogeneous  dielectrics.  By  using  a  s ‘ratified-medium 
foi  mulation  we  showed  that  these  poles  of  the  spectral  reflection  coefficient  satisfy  the  same  equation  as  the 
guidance  condition  in  inhomogeneous  waveguides.  Therefore,  in  terms  of  wave  numbers,  the  poles  are  the  same 
as  the  discrete  modes  in  the  waveguide.  The  radiation  modes  have  continuous  real  values  of  transverse  wave 
numbers  and  are  represented  by  the  branch  cut  on  the  complex  plane.  Based  on  these  results,  applications  of 
the  Gel’fand-Levitan-Marchenko  theory  to  optical  waveguide  synthesis  with  the  rational  function  representa¬ 
tion  of  the  transverse  reflection  coefficient  are  discussed. 


1.  INTRODUCTION 

The  electromagnetic  inverse  problem  is  to  find  unknown 
parameters  of  an  object  from  its  responses  to  given  elec¬ 
tromagnetic  signals.  The  electromagnetic  inversion 
method  is  an  emerging  technique  in  many  fields,  such  as 
geophysical  media  and  ionosphere  probing,  medical  imag¬ 
ing,  nondestructive  material  characterization,  and  remote 
sensing.  These  applications  belong  to  the  identification 
problem,  which  is  to  determine  the  unknown  parameters 
from  measured  data.  Recently,  inversion  theory  has  been 
applied  to  synthesis  (or  design)  problems  in  which  the  re¬ 
sponse  of  the  object  is  not  measured  but  prescribed  from 
design  criteria.  One  important  application  is  the  design 
of  dielectric  optical  waveguides.'"^  This  design  problem 
is  to  determine  the  refractive-index  profile  of  a  graded- 
refractive-index  guided  wave  device  for  given  require¬ 
ments  of  modal  structure.  One  of  the  advantages  of  using 
a  graded-refractive-index  optical  waveguide  is  that  the 
core  region  can  be  made  wider  relative  to  the  homoge¬ 
neous  guide.  Therefore  it  is  easier  to  fabricate  the  guide 
at  optical  wavelengths. 

For  a  one-dimensional  planar  problem,  as  shown  in 
Fig.  1,  the  dielectric  profile  e(x)  can  be  obtained  from  the 
reflection  coefficients  r  measured  at  the  surface  x  =  0. 
This  is  a  well-known  inverse-scattering  problem.  In  the 
spectral  inverse-scattering  theory  [Gel’fand-Levitan- 
Marchenko  (G-L-M)],  closed-form  solutions  of  «(x)  may 
be  obtained  if  the  reflection  coefficient  r  is  represented  by 
a  rational  function.'*'’  In  the  design  problem,  for  given 
requirements  of  possible  guided  modes  in  the  waveguide, 
the  first  step  is  to  determine  the  transverse  reflection  co¬ 
efficient  on  which  the  spectral  inverse-scattering  theory 
(G-L-M)  can  be  applied.  This  reflection  coefficient 
should  contain  all  the  modal  information  required  in  the 
slab.  The  second  step  is  to  apply  the  inverse  theory  to 
solve  the  index  profile.  This  step  has  been  discussed  at 
great  length  in  the  literature,  while  the  first  step  has  re¬ 
ceived  little  discussion.  An  alternative  design  procedure 
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is  to  solve  the  index  profile  directly  from  the  modal  re¬ 
quirements  without  constructing  a  reflection  coefficient. 
That  is  not  discussed  in  this  paper.  The  motivation  for 
using  a  reflection  coefficient  here  is  to  utilize  the  previ¬ 
ous  results  in  inverse  scattering,  especially  the  rational 
function  representation  of  the  reflection  coefficient  for 
closed-form  solutions.  The  purpose  of  this  paper  is  to  in¬ 
vestigate  the  feasibility  of  designing  optical  waveguides  by 
constructing  a  transverse  reflection  coefficient  as  a  start¬ 
ing  point. 

There  are  two  types  of  measurement  in  one-dimensional 
frequency-domain  problems,  frequency  diversity  and  angle 
diversity.  These  correspond  to  the  two  independent  vari¬ 
ables,  the  angular  frequency  w  and  a  component  of  the 
wave  vector.  If  a  layered  model  is  used  to  describe  the 
inhomogeneous  medium,  owing  to  the  phase-matching 
condition  on  the  boundaries  between  layers,  k.  will  be  in¬ 
dependent  of  X,  i.e.,  all  the  layers  have  the  same  k,.  Since 
ky  =  0,  there  are  only  two  components  in  the  wave  vector 
k,  and  ky.  The  components  of  a  wave  vector  satisfy  the 
dispersion  relation  in  each  layer.  Therefore  only  one 
component  is  independjent.  In  region  €i  in  Fig.  1  we 
choose  ku  as  the  independent  variable.  The  transverse 
reflection  coefficient  is  written  as  a  function  of  amd  to, 
i.e.,  r{ki„to),  which  is  consistent  with  the  notation  in  pre¬ 
vious  work,'  where  r(A,,)  is  written  as  r(k).  In  optical- 
waveguide-design  problems,  since  the  operating  frequency 
is  usually  fixed,  riku,  to)  reduces  to  r(ku).  This  is  differ¬ 
ent  from  the  profile  identification  problem  in  Refs.  5  and 
7,  where,  for  fixed  normal  incidence,  ku  =  k  =  to/c,  rico)  is 
actually  considered.  In  this  paper  we  discuss  only  r(ku^ 
for  the  waveguide-design  problem.  Furthermore,  the  ra¬ 
tional  function  representation  of  r{ku\  where  the  (K)les 
of  the  reflection  coefficient  are  related  to  the  electro¬ 
magnetic  modes  supported  by  the  waveguide,  is  consid¬ 
ered.  The  correspondence  between  the  poles  and  the 
modes  is  studied. 

An  integral  formulation  known  as  the  Sommerfeld  inte¬ 
gral"'"’  is  usually  used  in  calculating  the  fields  of  a  dipole 
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Fig.  1.  One-dimensional  inverse-scattering  problem. 

in  layered  media.  It  is  a  well-known  result  that  the  poles 
of  the  Sommerfeld  integrand  correspond  to  the  modes  in 
the  media.  This  result  can  be  obtained  from  the  residue 
theorem  in  the  theory  of  complex  variables.  For  instance, 
the  Green’s  function  in  a  layered  medium  can  be  writ¬ 
ten  as 

G(r,  r')  =  jj  dks[A(ks)  +  r(k,)B(kJ] 

X  exp[(ks  •  (Fj  -  r/)],  (1) 

where  =  yy  +  zz,  dk,  =  dk^dk;,  and  k,  =  yi,  +  zk^  is 
the  wave  vector  in  the  plane  of  the  layers  (Fig.  1).  k,  re¬ 
duces  to  zk.  since  ky  =  0.  A(k,,)  and  BCk,)  are  functions 
related  to  the  medium.  Surface-wave  modes  are  the  poles 
of  the  Green’s  function,  which  are  from  the  poles  of  the 
transverse  reflection  coefficient  here  is  chosen  as 
the  independent  wave  number.  Studying  the  Green’s 
function  in  the  Sommerfeld  integral  is  one  way  to  obtain 
the  relation  between  the  poles  of  r  and  the  modes  in  the 
medium.  In  constructing  a  rational  reflection  coefficient 
in  the  design  problem,  it  is  critical  to  relate  the  guidance 
condition  to  the  pole  condition  of  the  reflection  coeffi¬ 
cient.  To  understand  the  pole  condition  of  r  in  terms  of 
the  field  distribution,  in  this  paper  we  propose  another  ap¬ 
proach,  which  is  based  on  the  analysis  of  the  electromag¬ 
netic  reflection  coefficient  and  the  guidance  condition  of 
a  layered  medium. 

For  a  given  inhomogeneous  medium,  the  reflection  coef¬ 
ficient  is  derived  in  recursive  form  by  using  the  layered- 
medium  model.  The  pole  condition  for  the  reflection 
coefficient  can  then  be  obtained.  A  mode  is  a  possible 
field  structure  in  certain  physical  systems;  electromag¬ 
netic  modes  are  the  possible  solutions  of  Maxwell’s  equa¬ 
tions  with  certain  boundary  conditions.  In  a  layered 
medium  the  guidance  condition  is  a  necessary  condition 
for  the  existence  of  guided  modes  derived  from  the  bound¬ 
ary  conditions.  We  show  that  the  pole  equation  for  the 
reflection  coefficient  is  the  same  as  the  guidance  condi¬ 
tion  in  the  medium.  Thus  the  poles  of  the  reflection  coef¬ 
ficient  are  indeed  the  modes  of  the  guiding  structure. 
The  starting  point  of  this  study  is  a  single-layer  medium 
that  is  a  homogeneous  dielectric  slab  waveguide.  In  this 
simple  case  the  possible  solutions  for  the  pole  equation  can 
be  fully  analyzed.  For  a  general  iV-Iayer  medium,  a  proof 
is  given  that  the  guidance  conditions  in  all  the  layers  are 
equivalent.  Furthermore,  this  guidance  condition  is 
.■"hown  to  be  equivalent  to  the  pole  equation  of  the  total 


reflection  coefficient.  The  solutions  of  the  guidance  con¬ 
dition  give  both  the  guided  and  the  leaky  modes.  Radia¬ 
tion  modes  are  not  poles  of  the  reflection  coefficient.  The 
results  on  the  pole-mode  relation  are  used  to  interpret  the 
poles  of  the  rational  reflection  coefficient  in  the  inverse¬ 
scattering  problem  applied  to  optical  waveguide  design. 


2.  HOMOGENEOUS  DIELECTRIC 
WAVEGUIDE 

To  find  the  relationship  between  the  modal  structure  and 
the  transverse  reflection  coefficient,  we  first  analyze  the 
homogeneous  dielectric  waveguide  shown  in  Fig.  2.  Al¬ 
though  this  is  the  simplest  case  of  general  dielectric  wave¬ 
guides  and  has  been  discussed  extensively,  it  serves  the 
purpose  of  illustrating  the  results  in  an  analytic  form. 

The  propagation  direction  is  -f  z  in  Fig.  2.  The  x  direc¬ 
tion  is  called  the  transverse  direction.  The  reflection 
coefficient  r(^i„a»)  is  obtained  by  solving  the  airect  scat¬ 
tering  problem  for  an  incident  plane  wave,® 


R-,0  +  Rni  exp(i2k,d) 
1  R-iojRoi  exp(i2k,d) 


(2) 


where  k,  is  the  wave  number  along  the  z  direction,  w  is  the 
angular  frequency,  d  is  the  thickness  of  the  slab,  and  R,j  is 
the  Fresnel  reflection  coefficient  for  a  plane  wave  incident 
from  medium  ii)  to  medium  (J).  Since 


i?-io  ~  ~  Ro 


(3) 


the  poles  are  determined  by 

7?oi^  exp(i2M)  =  1  (4) 

for  the  transverse-electric  (TE)  wave  case,  where 


Roi 


kx  -  ku 
kx  +  kix 


(5) 


The  dispersion  relations 

A  J  +  ku^  =  (6) 

where  feu  =  fei.  and 


k/  +  fe,2  =  fe2  (7) 

need  to  be  used  in  solving  Eq.  (4);  k,  is  the  same  in  all 
media  owing  to  the  phase-matching  condition.  It  is  easy 


Fig.  2.  Mode  diagram  for  discrete  modes  in  slab  waveguide. 
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to  prove  that 

r(/5:u  =  0)  =  -1,  (8) 

which  is  the  normalization  condition  of  the  rational  re¬ 
flection  coefficient. 

The  poles  of  are  those  values  of  that  satisfy 

Eq.  (4).  At  the  poles  r  ^  which  indicates  a  self¬ 
resonant  structure.  When  a  plane  wave  is  incident  upon 
the  slab  with  real  kix,  owing  to  conservation  of  energy 
|r|  s  1.  This  is  because,  for  real  values  of  and  kx, 
(/?oil  —  1,  the  pole  condition  will  not  be  satisfied.  The 
poles  of  rikix)  are  complex. 

In  Fig.  2  the  guided  modes  are  defined  as  the  waves 
propagating  along  the  +z  direction  and  exponentially  de¬ 
caying  along  the  transverse  direction  outside  the  slab. 
The  modes  in  the  dielectric  waveguide  are  determined 
from  the  boundary  conditions.  For  a  TE  wave,  the  E, 
fields  in  the  three  regions  can  be  expressed  as 

E-i  =  E-i  ex’p(ikxz)exp{.-ik .xx) 

for  region  (-l),x  <  0,  (9a) 
iJo  =  exp(iA!j2)[A  expdAjjc)  -t-  B  exp(-z7!,x)] 

for  region  (0),0  <  x  <  d,  (9b) 

El  =  El  exp(.ikxZ)exp(ikixX)  for  region  (l),x  >  d. 

(9c) 

At  X  =  0  in  region  (0),  the  right-going  wave  is  due  to  the 
reflection  of  the  left-going  wave  at  the  boundary,  which 
yields  A  =  BRo,-u-  With  floi-n  =  -R-io  and  using 
Eq.  (3),  we  have 

A  =  BRoi .  (10) 

At  X  =  (f  in  region  (0),  the  left-going  wave  is  due  to  the 
reflection  of  the  right-going  wave  at  this  boundary;  thus 

B  exp{-ik,d)  =  Aexp{ikxd)Riji,  (11) 

and  combining  Eqs.  (10)  and  (11)  yields 

exp{2ikxd)Roi^  =  1.  (12) 

Equation  (12)  is  a  transverse  resonance  condition  satisfied 
by  any  field  distribution  of  Fig.  2  (as  we  see  in  matching 
boundary  conditions)  and  is  known  as  the  guidance  condi¬ 
tion.*'”  This  equation  is  the  same  as  Eq.  (4).  Thus  the 
guidance  condition  is  equivalent  to  the  pole  condition  for 
the  reflection  coefficient.  Equation  (12)  can  also  be  ob¬ 
tained  by  using  the  boundary  conditions  at  x  =  0  and 
X  =  rf  for  both  the  E  and  the  H  fields  from  Eq.  (9). 

For  a  fixed  real  angular  frequency  co  there  is  only  one 
unknown  in  Eq.  (12).  This  unknown  can  either  be  kix  or 
kx-  The  possible  solution.s  that  satisfy  the  guidance  condi¬ 
tion  [Eq.  (12)]  can  be  discussed  in  the  following  three 
cases  of  k,: 

Case  1.  kx  is  Real 

(a)  kix  is  also  real.  This  gives  lexp(i2A,<i)l  =  1  and 
<  1.  This  case  cannot  satisfy  guidance-condition 
equation  (12).  We  assume  that  e  >  e,.  For  cj  >  e,  since 
kix^  =  kx^  +  ki^  -  k^,  k,  and  ^1,  will  both  be  real,  as  in  the 
present  case.  Thus  there  is  no  discrete  mode  for  ti  >  e. 
(b)  kix  is  purely  imaginary.  Here  k\^  -  k'^  =  ki^  ~ 


<  0,  Cl  <  €,  and  kix  =  kx  +  ki  -  is  real.  Thus,  if 
kix  is  complex,  ku  has  to  be  purely  imaginary.  In  this 
case  |exp(z2^,cf)|  =  1,  |i?oil  =  1,  so  that  it  is  possible  for 
the  guidance  condition  [Eq.  (12)]  to  have  a  finite  number 
of  solutions  with  discrete  kx  values,  kx^  =  ku^  +  k^  - 
ki^  <  k^  —  ki^.  The  proper  kix  values  are  on  the  positive 
imaginary  axis  on  the  complex  kix  plane.  Since  k^  is 
positive  imaginary,  the  wave  is  decaying  along  the  x  direc¬ 
tion  outside  the  slab.  There  are  improper  A:  i,  solutions  on 
the  negative  imaginary  axis  of  the  ku  plane.  The  solution 
of  Eq.  (12)  in  the  present  case  also  includes  the  trivial 
root  kx  =  0,  which  gives  a  zero  field  everywhere.  From 
Eq.  (10),  if  kx  =  0,  /?oi  =  “1  and  then  A  =  -B.  From 
Eq.  (9),  Eo  =  0.  Since  E,  is  continuous  at  the  boundaries, 
E-i  =  0  and  Ei  =  0. 


Case  2.  k,  is  Purely  Imaginary 
This  implies  that  k^^  >  k^  and  k^  is  real. 

(a)  e  >  ei,k  >  ki,  so  that 

kix^=kx^-{k^-  ki^)<kx^<0,  (13) 

and  ku  is  also  an  imaginary  number.  Let  A,  =  ia  and 
ku  =  ioi  for  decaying  waves  outside  the  slab,  and  a,  >  0 
and  |ai|  >  |a|. 

(i)  a  >  0,  |floil  <  1.  and  exp(-2ac?)  <  1,  which  is  not  a 
possible  solution  to  Eq.  (12). 

(ii)  a  <  0,  |i?oil  >  1.  and  exp(-2ad)  >  1,  which  is  also 
not  a  possible  solution  to  Eq.  (12). 

(b)  €  <  ei,k  <  ki. 

(i)  is  real,  ki  >  k,.  ii?oil  =  L  and  |exp(-2ac()|  >  1 
or  <  1.  This  is  not  a  possible  solution  for  Eq.  (12). 

(ii)  ku  is  imaginary,  ki  <  k,,  the  same  as  in  case  2(a). 

From  case  2  it  can  be  concluded  that  purely  imaginary 

kx  cannot  be  the  poles  of  the  reflection  coefficient. 


Case  3.  kx  Has  Both  Real  and  Imaginary  Parts 
Since  =  k'  -  kx^,  k^  is  also  complex,  and  kx  =  kx'  + 
ikx".  Since  the  propagation  direction  is  +z,  k,  >  0  and 
kx"  >  0.  Defining  kx  —  kx  +  ikx"  and  from  kx^  =  k^  - 
kx^  =  k'^  -  kx'^  +  k,"^  -  i2kx'kx”  we  get  kx'kx”  =  -kx'k." 
For  the  same  reason,  ku  =  ku'  +  iku",  and  ku^  =  ki^  - 

kx^  =  ki^  -  kx'^  +  kx"'  -  i2kx'kx"  and  ku'ku"  =  —kx'kx". 

Since  kx'kx"  >  0  and  kx'kx"  <  0,  ku'ku"  <  0.  Outgoing 
waves  require  that  ku'  >  0  in  regions  (-1)  and  (1);  hence 

ku"  <  0. 


There  are  two  equivalent  cases  for  >  OandA,'  <  0. 
We  choose  kx'  >  0  and  kx"  <  0.  This  means  that 
lexp(-2A,'y)|  >  1.  IfEq.  (12)  is  satisfied,  |i?oi|  <  1-  Since 


l^o.l 


ikx'  -  ku')’  +  (kx"  -  ku")^ 
(kx'  +  ku')^  +  ikx"  +  ku")' 


(14) 


gives  |/?oil  <  1,  this  is  a  possible  solution.  These  solutions 
are  called  the  leaky  modes  since  the  power  leaks  away 
from  the  surface.  As  is  illustrated  in  Fig.  3(b),  outside 
the  slab  the  wave  is  growing  away  from  the  waveguide 
along  the  x  direction  but  decaying  along  the  +z  direction. 
In  this  mode  diagram,  the  total  wave  is  propagating  along 
the  real  part  of  the  wave  vector,  k,/?,  and  decaying  along 
the  imaginary  part  of  the  wave  vector,  ku. 

In  conclusion,  the  possible  solutions  to  the  guidance 
condition  [Eq.  (12)]  are  discrete.  These  solutions  can  be 
classified  into  two  cases: 
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(a)  gutded  modes 


(c)  radiation  modes 

Fig.  3.  Wave  vectors  in  slab  waveguide  for  (a)  guided  modes, 
(b)  leaky  modes,  and  (c)  radiation  modes. 

(1)  Guided  modes,  k,  is  real,  and  is 

purely  imaginary  with  a  positive  imaginary  part. 

These  are  also  called  the  surface-wave  modes  because 
the  wave  is  guided  inside  the  two  surfaces.  Outside  the 
surfaces  the  wave  decays  exponentially.  This  surface 
wave  is  associated  with  the  total  internal  reflection  at  the 
surface  boundaries. 

(2)  Leaky  modes,  k,  is  complex,  and  ki,'  >  0, 

ki."  <  0. 

There  are  an  infinite  number  of  leaky  modes  that  are 
also  discrete.  They  usually  do  not  contribute  to  the  com¬ 
plete  eigenfunction  set  of  field  solutions  of  the  open 


waveguide.  Leaky  modes  lie  on  the  improper  Riemann 
sheet  [which  is  separated  from  the  proper  Riemann  sheet 
by  the  cut  of  Im(Au)  =  0]  of  the  complex  plane.  Most 
of  them  will  not  be  excited.  A  finite  number  of  leaky 
modes  may  contribute  to  the  total  field  when  the  excita¬ 
tion  condition  is  satisfied  (Ref.  8,  p.  326).  In  this  case  the 
leaky-mode  amplitudes  will  not  diverge. 

The  cutoff  frequency  of  a  guided  mode  is  determined  by 
the  condition  ^  -  ki^.  At  cutoff,  A,  = 

fei,  =  0,  and  Roi  =  1.  From  the  guidance  condition 
[Eq.  (12)],  k,  =  rntr/d  for  m  =  0,  ±1,  ±2, _  Hence 

,  mir  1 

“■  t  .  V 1  2  ’ 

d  (fr  -  eir) 

The  cutoff  frequencies  for  the  TE„  [also  true  for 
transverse-magnetic  (TM„)j  modes  are 

.  _  cm  1 

W  (€,  -  eix)*'' 

As  the  frequency  decreases,  when  f  <  /|.„,  the  mth  guided 
mode  will  be  cut  off.  As  a  solution  to  Eq.  (12)  this  mode 
becomes  a  leaky  mode  that  has  both  real  and  imaginary 
parts  in  k^.  The  leaky  mode  is  different  from  a  radia¬ 
tion  !.jode. 

The  radiation  modes  are  defined  as  possible  field  distri¬ 
butions  in  the  slab  with  real  ku  in  regions  (-1)  and  (1). 
These  modes  carry  energy  infinitely  far  away  and  are  not 
guided  locally.  In  terms  of  the  spectral  parameters  k^, 
the  difference  between  radiation  modes  and  leaky  modes 
is  that  leaky  modes  have  an  imaginary  part  kx^'  <  0  and 
radiation  modes  have  a  continuous  real  kxx- 

As  discussed  above,  it  is  impossible  to  have  fields  of 
Eqs.  (9)  with  real  ki^;  hence  the  radiation  modes  are  not 
possible  with  the  field  structures  of  Fig.  2.  The  radiation 
modes  are  not  poles  of  the  reflection  coefficient  r{kix,to). 

Another  possible  field  structure  is  shown  in  Fig.  4,  where 

£-1  =  E-i  exp(ikx^)exp(iklJ,x)  +  r(ki,)E.i  expiik^z) 

X  exp(-z/ti,x) 

for  region  (-1),  x  <  0,  (15a) 

£o  =  exp(!^xZ)[A  exp(iA:,x)  +  B  exp(-i7E,x)] 

for  region  (0),0  <  x  <  d ,  (15b) 

El  =  El  exp{ikxZ)exp(ikixX)  for  region  (l),x  >  d, 

(15c) 


Fig.  4.  Mode  diagram  for  radiation  modes  in  slab  waveguide. 
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Note  that  Re(Ai,)  >  0,  and  the  field  in  region  (-1)  has 
an  incident  term  that  is  different  from  the  one  in  Fig.  2. 
With  this  field  distribution,  at  x  =  d,  B  exp(—ikjd)  = 
A  e\p(ikjd)Ro\  is  still  valid.  However,  at  *  =  0,  A  = 
BRoi  +  TioE-i,  which  is  different  from  Eq.  (10).  This  can 
also  be  checked  by  matching  the  boundary  conditions  on  E 
and  H  at  x  =  0.  Equation  (12)  does  not  hold  in  this  case, 
so  that  ku  can  be  real  numbers;  therefore  the  radiation 
modes  exist.  The  modal  structure  of  Fig.  4  is  equivalent 
to  putting  a  source  in  region  (-1)  so  that  the  wave  is  inci¬ 
dent  from  X  =  -  00.  Now  the  k  i,  value  can  change  continu¬ 
ously.  Radiation  modes  form  the  continuous  spectrum  of 
eigenfunctions  of  a  dielectric  waveguide.  Continuous  ra¬ 
diation  modes  along  with  the  discrete  guided  modes,  in 
general,  form  a  complete  eigenfunction  set  in  which  the 
field  from  an  arbitrary  source  can  be  expanded.®  The 
continuous  spectrum  also  includes  the  radiation  modes 
that  are  due  to  incident  waves  from  the  right-hand  region 
(1).  The  analogy  in  quantum  mechanics  is  that  the  un¬ 
bound  states  are  relevant  to  scattering  problems;  such 
problems  characteristically  involve  a  wave  that  is  incident 
from  the  exterior  region.”  Here  -  A,®  >  0, 

and  there  are  two  cases  to  be  considered; 


3.  GENERAL  INHOMOGENEOUS 
WAVEGUIDE 

In  this  section  we  analyze  a  general  inhomogeneous  wave¬ 
guide,  using  a  layered  model.  The  results  of  Section  2  are 
generalized  from  the  one-layer  slab  to  an  Af-layer  medium 
shown  in  Fig.  6.  Assuming  that  a  TE  wave  is  incident  in 
region  (-1)  from  the  left,  the  reflection  coefficient  can  be 
obtained  in  a  recursive  form  (Ref.  8,  p.  130,  Sec.  3.3;  note 
that  here  we  use  the  coordinate  system  as  in  Sec.  3.5)  as 

r{ku,(D)=^^<  (16) 

A-i 


where  A-i  and  B-i  are  the  wave  amplitudes  shown  in 
Fig.  6: 


r(ku,oj)  =  - —  -t- 
ti-io 


(1  -  1/R-io®) 
l/R_io  +  Bo/Ao 


(17) 


The  poles  are  the  zeros  of  the  pole  equation 


Cased)  is  real.  <  Ai,  and  A:.,  <  Ai. 

This  corresponds  to  the  case  in  which  the  power  is  inci¬ 
dent  from  the  exterior  regions,  (-1)  or  (1)  in  Fig.  4. 

Case  (2)  k,  is  purely  imaginary.  This  means  that 

^  Ail. 

This  case  is  sometimes  referred  as  the  evanescent 
mode^®®  since  the  wave  is  exponentially  decaying  along 
the  propagation  direction  +z. 


(18) 


Matching  the  boundary  conditions  of  E  and  H  at  x  =  <ii 
gives 

Bo  exp(i2k4\)  ^  (1  -  l/i?oi®)exp[t2(Aiu  +  Aijdi] 

Ao  Boi  1/Boi  exp(i2Aurfi)  +  Bi/Ai 


The  radiation  modes  lie  on  the  Sommerfeld  branch  cut. 
The  same  results  are  derived  from  the  Green’s-function 
approach.  These  radiation  modes  are  illustrated  in 
Fig.  3(c).  Note  that  k;  is  always  orthogonal  to  k/j,  which 
can  be  proved  by  calculating  k;  •  k^  and  using  k/k"  == 
so  that  k/  •  k/i  =  0. 

In  summary,  the  modal  structure  of  a  dielectric  slab 
(Fig.  5)  shows  the  pole  positions  in  the  complex  plane  cor¬ 
responding  to  the  modes  for  a  fixed  real  In  Fig.  5,  only 
half  of  the  complex  planes  are  shown,  owing  to  the  defini¬ 
tion  of  Re(^,)  s  0,  Re(Aj)  ^  0,  and  Re(Au)  s  0.  Note  that 
the  radiation  modes  lie  on  the  Sommerfeld  branch  cut  and 
that  the  leaky  modes  are  on  the  improper  Riemann  sheet. 

For  numerical  examples,  the  guided  and  the  leaky  modes 
are  calculated  in  the  slab  waveguide  for  two  different  fre¬ 
quencies.  The  first  one  is  at  microwave  frequency  f  = 
12  GHz,  d  =  \  cm,  =  10,  and  =  1.  Here  ko  = 
251  m“'  and  (A:®  -  Aii®)'  ®  =  754  m"'.  For  the  TE  polariza¬ 
tion,  by  calculating  the  cutoff  frequencies,  =  5.0m  GHz, 
we  find  three  guided  modes.  These  are  TE,>,  TEi,  and 
TE2.  Table  1  shows  the  calculated  wave  numbers,  Aj„  Aj,^, 
and  k^,  for  the  guided  modes  and  some  leaky  modes. 

Another  example  is  for  the  optical  frequency,  f  = 
2.255639  X  10”  Hz,  corresponding  to  Ao  =  1.33  /iim.  The 
parameters  are  d  =  2.5  p.m,  tr  =  12.96  (n  =  3.6), 
e,.  =  12.6025  (n,  =  3.55),  ko  =  4.7242  X  10®  m'*, 
(k‘‘  -  /fe,®)'®  =  2.8267  X  10®  m  ',  ki  =  16.78  X  10®  m’', 
and  k  =  17.02  X  10®  m"'.  The  cutoff  frequencies  are 
frm  —  1.0  X  lO^m  Hz.  The  results  are  shown  in  Table  2. 


Here  the  poles  of  Bo/Ao  are  not  the  poles  of  r(ku,  w)  be- 


0:  guided  mode 
x:  leaky  mode 


(b)  (c) 

Fig.  5.  Discrete  and  continuous  modes  shown  in  complex  (a)  k„ 
(b)  and  (c)  Ai,  planes. 
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Table  1.  Modes  of  Slab  Waveguide  at  12  GHz 


Modes 

k,{ia  ') 

Aii(m'‘) 

Guided 

TEo 

247.36 

(712.81 

755.87 

TE, 

487.74 

(575.66 

628.20 

TE, 

702.79 

(274.51 

372.30 

Leaky  (examples) 

1 

892.53  -  (130.58 

511.96  -  (227.64 

253.53  +  (459.69 

2 

1530.85  -  (271.30 

1340.39  -  (309.85 

315.16  +  (1317.84 

3 

3745.05  -  (458.58 

3669.45  -  (468.02 

469.11  +  (3660.96 

Table  2.  Modes  of  Optical  Slab  Waveguide 


Modes 

*.(10'  m'‘) 

AudO®  m-') 

A,(10®  m*') 

Guided 

TEo 

0.9749 

(2.6532 

16.9909 

TE, 

1.9170 

(2.0772 

16.9106 

TE, 

2.7267 

(0.7450 

16.7990 

Leaky  (examples) 

1 

3.5616  -  (0.6090 

2.2858  -  (0.9489 

16.6537  4-  (0.1302 

2 

4.8427  -  (0.9323 

3.9845  -  (1.1331 

16.3443  +  (0.2762 

3 

6.1177  -  (1.1450 

5.4565  -  (1.2837 

15.9286  +  (0.4398 

4 

7.3889  -  (1.3077 

6.8474  -  (1.4111 

15.3997  +  (0.6274 

5 

8.6573  -  (1.4406 

8.1975  -  (1.5214 

14.7474  +  (0.8457 

cause  at  these  values  r  is  finite.  For  an  arbitrary  layer  I, 

^  _  exp(t2fe<,df4i) 

■Al  Rui-fit 

^  [1  -  l/R;n^ii^3exp{t2[A,i<.n,  +  k,.]di..i} 
exp[i2*,/.i„d,*i]  + 

Al  _  exp(.-i2kiidi) 

Bi  Rm-v 

[1  -  l/Ri(;-ii^]exp{-t2[A,<-ii,  +  ki,]di) 

1/Riii-v  exp[-i2kii-i,^di]  +  Ai^JBi-i 
For  the  last  layer  n  =  <  —  1, 

7^  =  Rn,  exp{i2kn.d„) . 

An 

A  direct  result  from  Eq.  (17)  with  R-io  =  -1  is  the  nor¬ 
malization  condition  for  general  inhomogeneous  media, 

r(Au  =  0)  =  -l.  (23) 

The  modes  in  the  medium  are  studied  by  assuming  a 
TE  wave  with  A-i  =  0  in  Fig.  6.  A  guided  wave  has  a  real 
k,  and  decays  exponentially  in  regions  (0)  and  (t)  away 
from  the  waveguide;  its  energy  is  confined  in  the  layered 
region,  with  the  thickness  of  the  inhomogeneous  region  d, 
assumed  to  be  finite.  If  the  e(x)  approaches  a  constant  as 
X  -*  ae,  an  infinitely  extended  medium  can  be  reduced  to 
an  equivalent  finite-thickness  medium. 


(20) 


(21) 


(22) 


The  E,  fields  in  different  regions  can  be  written  as 
E-i  =  B-i  exp[ik,z)exp(-ik.ux) 

for  region  (-1), x  <  0,  (24a) 

Eg  =  exp(iA,2)[Ao  exp(ik,x)  +  Bo  exp{-ik,x) 

for  region  (0),0  <  x  <  <fi,  (24b) 
El  =  exp{ikiZ)[Ai  expdki^x)  +  Bj  exp(-tAux) 

for  region  (l),di  <  x  <  d2,  (24c) 

E„  =  exp(iM)[A„  expdknjx)  +  B„  exp(-iA„x) 

for  region  (n),  dn  <  x  <  d,,  (24d) 
E,  =  A,  exp(i^,z)exp(jA,,x) 

for  region  (t),  x  >  d,.  (24e) 

Z  wave 


Fig.  6.  Mode  diagram  in  inhomogeneous  medium. 
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The  necessary  condition  for  the  existence  of  any  field 
distribution  in  Fig.  6  is  the  guidance  condition  in  each 
layer.  For  layer  (/), 


RJRJ  =  1,  (25) 

where  Rj  =  BJAi  and  RJ  =  Ai/Bi  are  in  recursive  forms 
[Eqs.  (20)  and  (21)]. 

For  /  =  0,  the  guidance  condition  is 

R.°RJ'=1.  (26) 


Matching  the  boundary  conditions  for  the  E  and  H  at  the 
boundaries  gives 


Bo  _  exp(i2Aidi) 

Aa  Ra\ 

^  (1  -  l//?ni^)exp[t2(fe|,  +  fe,)c(i] 
1/Roi  exp(i2fti,di)  +  Bi/Ai 


Bo 


=  -R 


-10  • 


(27) 

(28) 


It  is  easy  to  see  that  the  guidance  condition  for  layer  (0)  is 
the  same  as  the  pole  equations  (18)  and  (19). 

Next  we  show  that  all  the  other  layers  have  the  same 
guidance  condition  or,  equivalently,  that  for  any  layers  I 
and  Z  -  1  the  guidance  conditions  are  the  same.  The 
guidance  condition  for  layer  (Z)  is 

Bi  ( exp{-i2ki,di) 

\  Rlil-\' 

+  [1  -  l/Z?;./-|.-]exp{-t2[fe,,-i'.t  +  A!,,]cZ;}\  ^ 
l/B,,,-„exp[-i2fe„.„.cZ,]  +  Ai-i/fii-J 

and  the  guidance  condition  for  layer  (Z  -  1)  is 

(exp[;2Z;,;-i.,cZ;] 

R.i-v, 

^  [1  -  l/Z?,,-i/]exp{t2[^;.  +  A;,1 

l//?,,-i,/ exp(i2^,,cZ/)  +  Bj/A/  /B(.i 

Although  Eq.  (30)  appears  different  from  the  condition 
in  layer  (Z),  Eq.  (29),  it  can  be  shown  that  Eqs.  (29)  and 
(30)  are  equivalent.  They  both  can  be  reduced  to  an  iden¬ 
tical  equation.  Equation  (30)  can  be  rewritten  as 

Bi/[1  -  l/B,/-i.i^]exp{Z2[^,^  + 

Ai\ -exp[i2A:,,-,,<Z,]/B.,-i,,  -)-  B,-JAi.x 

_ 1 _ V'  _ 

R,i-u  exp(i2ki,di)) 
which  can  be  further  simplified  as 


that  the  poles  represent  the  modes  in  the  inhomogeneous 
core  region  of  the  waveguide. 

Equations  (29)  and  (30)  are  the  conditions  that  any  field 
distribution  of  Fig.  6  must  satisfy.  Possible  solutions  of 
Eq.  (18)  can  be  studied  in  detail  as  in  the  single-slab  case. 
The  results  can  be  summarized  into  two  cases:  (1)  guided 
modes,  with  purely  imaginary  k^u  and  k,j,  and  (2)  leaky 
modes,  with  k-u  and  kt,  having  both  real  and  imaginary 
parts.  The  guided  and  the  leaky  modes  are  discrete. 
Again,  the  radiation  modes  do  not  satisfy  the  pole  equa¬ 
tion  [Eq.  (18)],  so  that  radiation  modes  are  not  poles  of  the 
reflection  coefficient  r(ki„(o).  Radiation  modes  require 
that  there  be  a  wave  incident  from  region  (  —  1)  or  region 
(t),  where  k-i,  and  are  real  numbers  that  can  vary  con¬ 
tinuously  from  0  to  Zt-i  and  0  to  k,,  respectively.  Since  the 
fields  in  one  layer  are  dependent  on  the  fields  in  all  the 
other  layers,  the  waves  are  considered  to  be  guided  in  all 
the  layers  rather  than  in  a  single  layer.  For  a  continu¬ 
ously  varying  eix),  the  layered  model  can  also  be  used  by 
making  each  layer  sufficiently  thin.  All  the  above  results 
are  valid  for  continuous  inhomogeneous  media. 


4.  IMPLICATIONS  FOR  INVERSE  PROBLEMS 

Inverse-scattering  problems  are  concerned  with  the  recon¬ 
struction  of  the  physical  properties  of  unknown  objects 
from  information  contained  in  their  scattering  data.  In 
contrast,  direct-scattering  problems  are  concerned  with 
determining  the  scattered  fields  from  known  physical 
properties.  In  inverse  problems,  the  scattering  data 
r{ku,  a>)  are  basically  a  set  of  measurements  that  relate  r  to 
ku  or  (o.  To  apply  the  G-L-M  inverse-scattering  theory, 
we  approximate  r(ku,u})  as  a  closed-form  function.  One 
example  is  that  the  reflection  coefficient  r{ki,,ui)  is  repre¬ 
sented  by  a  rational  function  of  ku  or  w,  and  r(Au,  at)  is 
prescribed  to  have  a  finite  number  of  poles  and  zeros; 
typically  for  a  good  approximation  many  poles  and  zeros 
may  be  needed,  similar  to  circuit  network  synthesis.'^ 
Closed-form  solutions  of  the  G-L-M  integral  equation  of 
inverse-scattering  theory  can  be  obtained  for  some  rational 
functions  of  r(a>). '  These  closed-form  solutions  are  impor¬ 
tant  not  only  in  providing  a  benchmark  for  other  numeri¬ 
cal  methods  but  in  solving  synthesis  problems. 

As  was  discussed  in  Section  1,  there  are  two  kinds 
of  inverse  profile  problem,  the  identification  and  the  syn¬ 
thesis  problems.  In  practical  profile  identification  prob¬ 
lems  it  may  be  difficult  to  determine  from  the  measured 
data  the  functional  form  of  r(k,j,w)  that  is  determined 
by  the  unknown  permittivity  profile.  Although  the  poles 
and  the  zeros  of  r(k\„aj)  can  be  obtained  by  using  Bode 
diagram  techniques  together  with  analytic  function  ex¬ 
tension  methods,  the  rational  function  approach  is  more 


B,  /[Bi-t/Ai-t)Ri.i.t,  exp{-i2[ki,  -)-  ^,,-i,JcZ;}  +  exp(-t2^/,<f,) 
A/\  Bi-JAi-i  exp[-i2k.i-i,.di  +  Ri.i-u] 


=  1. 


(32) 


Through  algebraic  simplification,  Eq.  (29)  can  also  be 
written  in  the  same  form  as  Eq.  (32).  Thus  Eqs.  (29)  and 
(30)  are  equivalent.  We  have  just  shown  that  for  all  layers 
Z  the  guidance  conditions  are  the  same.  This  guidance 
condition  is  the  same  as  the  pole  equation  [Eq.  (18)],  so 


appropriate  for  synthesis  problems.  The  synthesis  prob¬ 
lem  is  to  determine  the  permittivity  profile  needed  to 
produce  a  prescribed  reflection  coefficient  r{ki,,w).  In 
the  optical- waveguide-synthesis  problem,  the  relevant  pa¬ 
rameters  are  modeled  by  continuous  functions'^;  in  the 
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present  example  we  use  the  transverse  reflection  coeffi¬ 
cient  Since  r(ku)  will  be  in  general  combinations 

of  transcendental  functions,  which  is  evident  from  our 
discussion  of  direct  scattering,  we  rely  on  standard  elec¬ 
trical  engineering  practice'^  and  approximate  r(kix)  by  a 
rational  function  of  the  complex  variables  This  func¬ 
tion  can  be  specified  by  its  pole-zero  configuration.  The 
synthesis  procedure  is  summarized  as  follows:  (1)  The 
possible  modes  that  are  needed  in  the  waveguide  are 
specified.  These  determine  the  guided  wave  poles  of  the 
transverse  reflection  coefficient  Hkh).  A  rational  func¬ 
tion  r(ki^)  is  constructed  from  the  pole-zero  information. 
(2)  The  G-L-M  integral  equation  is  solved  to  obtain  the 
permittivity  profile  e(x).  The  modal  information  is  used 
to  verify  the  answer. 

The  study  of  the  correspondence  between  the  poles  of 
the  transverse  reflection  coefficient  and  the  modes  in  an 
inhomogeneous  medium  can  be  applied  to  interpret  the 
existing  results  in  inverse  scattering.  Previous  results®"’ 
are  from  solving  the  inverse  problem  for  r(ui)  data,  which 
is  an  identification  problem.  In  Ref.  5  a  two-pole,  no-zero 
rational  function  is  used.  These  examples  can  illustrate 
the  analytical  technique  but  cannot  be  used  directly  for 
optical  waveguide  design,  since  r(^i,)  rather  than  r(<a) 
should  be  used  in  the  optical-waveguide-design  problems. 
The  wave  equation  in  a  dielectric  medium  (Fig.  1)  at  a 
certain  frequency  to  is 

kMku,x)  =  0  forxaO,  (33) 

r>X 

which  is 

o' fiu~E(ku,x)  =  -  fAx)]E(kix,x} 

dx^ 

for  X  so.  (34) 

With  q(x)  =  ^o^[ei  ~  er(x)]  the  wave  equation  can  be  writ¬ 
ten  in  the  same  form  as  the  one  for  rito]  data.  The  inverse¬ 
scattering  theory  of  G-L-M®’  can  then  be  applied  directly 
to  r{ku)  data. 

The  procedure  using  the  inverse  method  to  design  an 
inhomogeneous  optical  waveguide  can  be  discussed  fur¬ 
ther.  If  the  guided  modes  are  specified  in  the  waveguide 
and  if  the  width  d  and  the  operating  frequency  oi  are  pre¬ 
scribed,  the  problem  is  to  determine  f(x).  In  general,  this 
problem  can  be  solved  only  by  numerical  methods.  If  the 
closed-form  results  obtained  by  applying  the  G-L-M  theory 
are  to  be  used,  the  first  step  in  constructing  a  transverse 
reflection  coefficient  riku)  from  the  modal  requirements 
is  critical.  This  is  also  the  main  concern  in  this  paper. 
So  far  we  have  discussed  only  part  of  this  construction, 
i.e.,  determining  the  poles  from  the  discrete  modes.  How 
to  determine  the  zeros  has  not  been  discussed.  For  the 
homogeneous  slab  waveguide  in  Section  2,  the  zeros  of 
rlkul  are  found  from  Eq.  (2),  which  gives 

R-io  +  Roi  e\p(i2k,d)  =  0.  (35) 

With  Eq.  (3)  this  gives  exp(i2k,d)  =  1,  so  that  k.  =  mr/d 

for  n  =  0,  il,  ±2 .  There  are  an  infinite  number  of 

zeros  of  rikiA-  For  a  general  inhomogeneous  waveguide 
in  Section  3,  the  zeros  of  riktA  are  found  from  Eq.  (17), 
which  yields 

=  -R-.o  =  Roi.  (36) 

An 
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Again,  this  equation  has  an  infinite  number  of  solutions. 
Note  that  zero-condition  equation  (36)  can  be  obtained  by 
multiplying  the  right-hand  side  of  pole-condition  equa¬ 
tion  (18)  by  Roi^.  This  implies  that  the  zero  condition  and 
the  pole  condition  are  related. 

As  we  saw  in  Section  3,  r{ku)  in  Eq.  (17)  is  not  a  rational 
function.  In  other  words,  r(ki,)  is  not  exactly  a  rational 
function  with  a  finite  number  N  of  zeros  and  a  finite 
number  M  of  poles; 

riku)  =  ro-jf - >  (37) 

n  iku  -  k,p) 

I 

where  ro  is  a  constant,  k,/s  are  zeros,  k,p’s  are  poles,  and  11 
denotes  multiplication.  Since  an  exponential  function  is 
involved  in  both  the  numerator  and  the  denominator,  rikui 
has  an  infinite  number  of  poles  and  zeros.  If  a  poly¬ 
nomial  expansion  (e.g.,  a  Taylor  expansion)  is  used  in  both 
the  numerator  and  the  denominator,  the  number  of  terms 
used  in  the  expansion  determines  the  number  of  zeros  of 
the  expanding  function.  Thus  the  representation  with  a 
finite  number  of  poles  and  zeros  of  riku)  is  only  an  ap¬ 
proximation.  It  is  impossible  to  include  the  infinite  num¬ 
ber  of  leaky  poles  in  riki,). 

The  inverse-scattering  method  gives  a  profile  of  e(x)  ex¬ 
tending  from  X  =  0  to  X  — >  -t-x.  Finite  width  d  of  the 
waveguide  will  require  a  truncation  of  the  profile.”  This 
is  another  approximation. 

The  following  implications  can  be  drawn  from  the  above 
discussions.  Rational  function  representations  of  the 
transverse  reflection  coefficient  r{ki^,w)  are  used  to  ob¬ 
tain  closed-form  solutions  of  the  G-L-M  integral  equa¬ 
tion.  Since  it  may  be  difficult  to  extract  the  pole-zero 
configuration  of  the  measured  reflection  coefficient,  the 
rational  function  approach  is  more  appropriate  for  the 
synthesis  problem.  In  the  application  to  optical  waveguide 
design,  riku)  data  need  to  be  prescribed  for  a  required 
operating  frequency  and  certain  modal  information.  The 
possible  modes  in  the  waveguide  determine  the  poles  of 
r(ku)-  Rational  function  riku)  is  an  approximation  of  the 
true  reflection  coefficient,  and  the  corresponding  inverse 
problem  can  be  solved  in  a  closed  form. 

5.  SUMMARY 

The  feasibility  of  applying  the  inverse-scattering  method 
with  a  rational  reflection  coefficient  to  optical  waveguide 
design  has  been  studied.  The  pole-mode  relation  in  the 
rational  reflection  coefficient  was  derived  by  using  a 
stratified-medium  model.  It  was  shown  that  the  pole  con¬ 
dition  for  the  transverse  reflection  coefficient  is  the  same 
as  the  guidance  condition  for  the  modes  in  a  general  inho¬ 
mogeneous  waveguide.  All  the  poles  of  the  reflection  co¬ 
efficient  correspond  to  the  discrete  modes,  which  are  the 
guided  and  the  leaky  modes.  The  radiation  modes  are 
continuous;  they  do  not  appear  as  poles  in  the  reflection 
coefficient  but  are  represented  by  the  Sommerfeld  branch 
cut  on  the  complex  k.  plane.  In  the  application  of  the  ra¬ 
tional  function  approach  to  optical  waveguide  synthesis, 
poles  in  the  rational  reflection  coefficient  riku)  deter¬ 
mine  the  possible  guided  and  leaky  modes  in  the  wave¬ 
guide.  To  design  the  waveguide  accurately  requires  a  large 
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number  of  poles  and  zeros  in  the  rational  function  repre¬ 
sentation  of  the  transverse  reflection  coefficient. 
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